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ABSTRACT 

Discretized numerical simulations are a powerful tool for investigation of non- 
linear MHD turbulence in accretion disks. However, confidence in their quanti- 
tative predictions requires a demonstration that further refinement of the spatial 
gridscale would not result in any significant change. This has yet to be accom- 
plished, particularly for global disk simulations. In this paper, we combine data 
from previously published stratified shearing box simulations and new global 
disk simulations to calibrate several quantitative diagnostics by which one can 
estimate progress toward numerical convergence. Using these diagnostics, we 
find that the established criterion for an adequate numerical description of lin- 
ear growth of the magneto-rotational instability (the number of cells across a 
wavelength of the fastest-growing vertical wavenumber mode) can be extended 
to a criterion for adequate description of nonlinear MHD disk turbulence, but 
the standard required is more stringent. We also find that azimuthal resolution, 
which has not often been extensively examined in previous studies, can signif- 
icantly affect the evolution of the poloidal magnetic field. We further analyze 
the comparative resolution requirements of a small sample of initial magnetic 
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field geometries; not surprisingly, more complicated initial field geometries re- 
quire higher spatial resolution. Otherwise, they tend to evolve to qualitatively 
similar states if evolved for sufficient time. Applying our quantitative resolution 
criteria to a sample of previously published global simulations, we find that, with 
perhaps a single exception, they are significantly under-resolved, and therefore 
underestimate the magnetic turbulence and resulting stress levels throughout the 
accretion flow. 

Subject headings: Black holes - magnetohydrodynamics (MHD) - stars:accretion 
- methods: numerical 



1. Introduction 



Numerical simulation is an important tool for understanding the dynamics and evolu- 
tion of accretion disks. In the last decade, the increase in computational power has made it 
possible to carry out global disk simulations of increasing complexity, as the state-of-the-art 



2003 ) to full three-dimensi onal general relativistic physics (jPe Villiers &: Hawley 



2003 



2009 



has r isen from pseudo-Newtonian dynamics (iHawley &: Krolikll200ll: iMachida fc Matsumoto 



2003 



Anninos et al.ll2005[) . most recently including toy- model thermodynamics ( iNoble et al. 



Gammie et al. 



Shafee et al.l 120081 ) . These simulations already have several achi evements to their 



credi t . Th ey have shown that the magneto-rotational instability (MRI; iBalbus fc Hawley 



199ll . Il998l ) can be effective in providing the required internal stress. They have shown that 
disks rapidly evolve to a near-Keplerian distribution of angular momentum regardless of the 
initial angular momentum distribution. They have demonstrated two mechanisms by which 
a large-scale field mag netic field can be atta ched to a black hole, either inward transport of 
truly large-scale field ( iBeckwith et al.l 120091) or spontaneous, if temporary, inflation of fiel d 
loops within the accretion flow (IMcKinney &: Gammid |2004J : iDe Villiers et al.l |2003| . |2005| ) . 
If the black hole rotates, these large-scale fields can support Poynting fiux-dominated jets. 

Nonetheless, a number of questions remain open regarding the quantitative quality of 
the predictions these simulations can make for real black holes in Nature. Some of these 
questions have to do with the adequacy of the physical approximations made, most notably 
the still rather crude account of thermodynamics in even the best of them. Radiation forces, 
although likely very important in numerous black hole accretion contexts, likewise remain 
largely omitted in global simulations. Others have to do with physical choices that must be 
made in order to calculate anything, but about which we have no real knowledge, particu- 
larly the initial strength and structure of the magnetic field. Confidence in the saturation 
strength of the field is greatest when it is statistically time-steady, yet very different from 
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its initial state. It is therefore desirable to begin with a weak magnetic field and let the 
natural dynamics of the situation strengthen it. Although there have been discussions in 
the literature that a certain initial field geometries are "most natural" (e.g., multiple loops 
Shafee et al.ll2008l ). all such arguments appear to be based on little more than subjective 
aesthetic judgments. We might hope that these choices have little effect on the outcome, 
but that must be demonstrated; in the case of jet-launchi ng, we already know that magnetic 
field geometry plays a major role ( iBeckwith et al.ll2008al ). 



Other questions have to do with the specifics of the numerical techniques employed. 
Every algorithm has certain strengths and weaknesses, and their impact depends on the 
nature of the problem being simulated. Geometric symmetry conditions may be imposed 
(e.g., imposing azimuthal boundary conditions periodic across a wedge rather than around 
a full circle in order to economize on computer time) whose consequences cannot easily be 
determined a priori. Transients due to arbitrary initial conditions must be eliminated. And 
every discrete numerical scheme must be shown to have fine enough resolution, both spatially 
and temporally, to describe adequately the physical processes involved. This last point is a 
special concern for accretion simulations because MHD turbulence is essential to the story, 
but true microphysical dissipation operates on a scale so many orders of magnitude smaller 
than the disk scale that no conceivable simulation can hope to resolve it. It is these more 
technical questions that are the province of this paper. 

In this paper, we will by no means attempt a full comparison of a broad range of 
codes and algorithms. We will, however, present some comparative data in the shearing-box 
context in order to estimate the potential level of contrast in saturated stress and other 
properties among codes utilizing different algorithms run at similar resolutions. 

We will also pass over th e second question, that of overly-restrictive symmetry condi- 
tions. We will only note that ISchnittman et al.l ( l2006l ) showed that in a full-27r simulation 
there is substantial power in stress fluctuations on azimuthal scales as long as ~ 1 radian. 
Given that flnding, wedges of at least a quarter-circle are certainly required. 

The third question, how to eliminate transients, poses somewhat different issues, de- 
pending upon whether the contex t is stratifled shear ing-boxes or global disks. In the former 
case, it has long been established (IHawley et al.lll995l ) that the principal transients are erased 
i n ~ 10 orbits, but t he turbulence exhibits signiflcant long-term, chaotic fluctuation power 
(IWinters et al.ll2003l ). even over durations as long as the very longest simulations, > 500 
orbits. Nonetheless, despite the long-term variability, it is possible to deflne time-averages 
reasonably well. We will discuss quantitative definition of transient removal in global disk 
simulations more fully in § [51 Here it suffices to say that statistical stationarity in global 
simulations has a dual meaning. One sense is the same as for shearing-boxes: given a surface 
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density and an orbital shear, the system should achieve a statistical steady-state in the am- 
plitude of the MHD turbulence. The other sense is that the local surface density must have 
a well-defined mean value over some period of time. This is more difficult, and can never be 
achieved over the entire simulation volume. 

The fourth and last question, how to determine whether a given simulation has adequate 
resolution to return accurate quantitative results, will occupy most of our attention. One 
expects convergence of the numerical solution to an exact solution as the grid size goes to 
zero, i.e.. Ax — )■ 0. Since that limit is never attained, convergence generally refers to the 
observation that a given quantity approaches some fixed value as Ax is reduced. Different 
quantities, of course, have different convergence rates (in practice, the formal convergence 
rate, determined by the order of the scheme, need not be manifest at a given finite resolution), 
and individual quantities may themselves be subject to several different criteria. The issue 
of convergence is further complicated when evolving equations without explicit resistivity 
and viscosity. In that case, dissipation occurs at the grid scale, and the effective Reynolds 
number of the system becomes a function of resolution. 

As an illustration, consider the central quantity of accretion physics, the local stress. Its 
saturation depends on the interplay between several different processes. First, small magnetic 
fluctuations are amplifled by the MRI; the spatial resolution must therefore be great enough 
that the fastest-growing linear modes — for both poloidal and toroidal fluctuations — e-fold 
at the correct rate. Next, nonlinear couplings between different wave modes must also be 
described properly. Adequate resolution for this process depends upon the lengthscale of the 
shortest important modes, but it is hard to determine a priori what that lengthscale may 
be. Finally the turbulenc e must be dissipated at scales much less than the stirring scale. 



Lesur fc Longarettil (120101 ) argue that in feasible simulations, all resolvable wavelengths are 
coupled together by non-linear processes. They also argue, however, that in real disks the 
range of non-local coupling in wavevector space is much smaller than the dynamic range 
between the stirring scale and the physical dissipation scale, permitting the creation of a 
genuine inertial range. 

Initial conditions may also play a role in determining the necessary resolution. Because 
the fastest-growing MRI wavelength is proportional to the fleld strength parallel to the 
wavevector, the weaker the initial fleld, the flner the spatial resolution must be. In the non- 
linear stage when the magnetic fleld has a flxed large-scale element, even modest resolution 
suffices to describe that portion reasonably well, although a full description of the turbulence 
places stronger demands. On the other hand, when even the largest scales maintained by 
the physics are smaller, the minimum necessary resolution becomes correspondingly finer. 
Popular initial conditions for zero net-fiux magnetic field configurations illustrate this effect. 



- 5 - 



When the initial magnetic field is a single set of nested poloidal loops, many of these field 
line loops become attached to the black hole horizon, forming a radially large-scale magnetic 
structure linking the black hole and the disk. On the other hand, when the initial field is 
quadrupolar, with a pair of nested poloidal loops on either side of the midplane, these loops 
are free to shrink in size, demanding a much finer spatial grid if excessive resistive losses are 
to be avoided. 

For reasons of practicality, simulators tend to choose resolutions that are as fine as pos- 
sible while still consistent with their computing budget. It is therefore often not feasible 
to test directly for convergence by performing a new simulation with more gridzones. Al- 
though some sense of the rate of convergence can be obtained by computing models at lower 
resolution, that procedure would not reveal the failure of the highest resolution simulation 
to resolve some important feature — it is entirely possible that the reason the results do not 
change appreciably with changing resolution is that even the best-resolved case is incapable 
of sustaining an important physical effect. 

Our approach to these problems is to begin with shearing-box calculations, in which 
it is generically easier to reach high resolution than in global disk simulations. We will 
search for quantities related to, but different from, the stress and that scale with resolution 
in a way that can be calibrated. With these in hand, it becomes possible to gain some 
idea of how far along the path toward convergence a global simulation has reached. We 
will apply these measures both to a number of older simulations and to some new ones 
we have carried out especially for this purpose. Full general relativity is not necessary to 
achieve these goals because the saturation of MHD turbulence is an essentially Newtonian 
process. Consequently, it is sufficient to carry out simulations using the simpler and less 
computationally demanding Zeus algorithm and a pseudo-Newtonian potential. 



2. Numerical Setup and Initial Conditions 

In this study we explore some of the effects of resolution and initial conditions first by 
examining the results from local stratified shearing box simulations, and then by computing 
a series of global accretion simulations. The shearing box results are, in all but two cases, 
taken from the literature. Our sample features several resolutions, box sizes, and numerical 
algorithms. Some of the numerical aspects of those simulations are described in ^ for a 
detailed account the reader should consult the papers where these simulations were originally 
presented. 

Here we describe the numerical setup for the global accretion simulations carried out 
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specifically for this paper. As we will be focusing on the body of the disk itself rather than 
the interactions with the black hole or any jets or winds that form, it is sufficient to use non- 
relativistic MHD and to work in the pseudo- Newtonian potential to approximate the gravity 
of a black hole. We evolve the equations of Newtonian MHD in cylindrical coordinates 

§^ + v-(H = o (1) 

Q / tj2 \ / B \ 

+ (Pv ■ V)v = -V P + Q + — - pV$ + — ■ V B (2) 



dt ' ' V Stt 7 V47r 



^ + V-(pev) = -(P+Q)V-v (3) 



c?B 

— = V X (v X B) (4) 

where p is the mass density, e is the specific internal energy, v is the fiuid velocity, P is 
the pressure, $ is the gravitational potential, B is the magnetic field vector, and Q is an 



explicit artificial viscosity of the form described by IStone fc Norman! (119921 ). To model a 
black hole gravitational field we use the pseudo-Newtonian potential of Paczyhski & Wiita 
(1980) which is 

$ = , 5 

r — Tg 

where r is spherical radius, and Vg = 2GM/c^ is the "gravitational radius," akin to the 
black hole horizon. For this potential, the Keplerian specific angular momentum (i.e., that 
corresponding to a circular orbit) is 

hep = {GMryl^^—, (6) 

r Tg 

and the angular frequency Vt = l/R^. The orbital period at a radius r is Porb = 2ttQ~^ = 
"g)/r. The innermost stable circular orbit (ISCO) is located at Vms = 3rg. We use 
an adiabatic equation of state, P = pe(T — 1) = Kp^ , where P is the pressure, p is the mass 
density, e is the specific internal energy, if is a constant, and T = 5/3. Radiation transport 
and losses are omitted. Since there is no explicit resistivity or physical viscosity, the gas 
can heat only through adiabatic compression or by artificial viscosity which acts in shocks, 
including weak shocks found in the turbulence. 

The code used is a time-explicit Eulerian finite-differencing Zeus code for MHD (Stone 
& Norman 1992a; Stone & Norman 1992b; Hawley & Stone 1995)0 We set GM = c = 1, so 



^Here the term Zeus refers to the algorithm rather than a specific code implementation. There are several 
publicly available versions of Zeus as well as a large number of individually developed versions, such as was 
employed here. 
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that rg = 2M. Time and distance units are given in terms of the mass M. 

The initial conditions consist of an orbiting gas torus with an angular momentum dis- 
tribution parameter q = 1.65 {Q oc R~'^), and K = 0.0034. The pressure maximum radius is 
R = 35M and the inner torus edge is at 20M. The orbital period at the pressure maximum is 
1227M. Using the definition of scale height H' in terms of the density moment, H' / R = 0.1 
at the pressure maximum. Explicitly fitting a Gaussian function oc exp[~{z/Hy] to the 
vertical pressure distribution gives a value of H/R ^ 0.16. Note that H = \/2Hg, where 
Hg is more commonly regarded as the Gaussian scale height. In an isothermal thin disk, 
Hg = Cs/^l. In internal code units the initial mass of the torus is 6096 (assuming a full 
2n azimuthal domain; the actual computational domain runs from to 7r/2). Of this total 
mass, 19% lies inside of i? = 35M. 

We use two initial field configurations and two average field strengths. The first is the 
standard dipole loop, in which the vector potential is written in the form 

= C{p - Pent) (7) 

so that the field lines run along contours of constant density. Here C is a constant that sets 
the overall field strength and is set to zero wherever p < p^ut- This configuration is referred 
to as "one- loop . " We initialize a "two-loop" simulation using the vector potential function 



of IShafee et al.l (120081 ) . namely 

= [{p - PcutY'-'"] ' sin [Hr/S)/T] (8) 

where r is the spherical coordinate radius, pcut is set at 20% of the density maximum (thus 
confining the initial field to well within the edge of the initial torus), S = l.lrj„, where 
Tin = 20M is the initial inner edge of the torus, and T = 0.16 (see Fig. [1]). For either 
geometry, the initial field strength is normalized to a /3 value, either 100 or 1000, defined 
as the ratio of the total volume-integrated gas pressure to the volume-integrated magnetic 
pressure. 

Boundary conditions are periodic in and outflow along the z and R boundaries. One 
consequence of cylindrical coordinates and our use of an inner boundary radius located at 
Rin > ^g, is that there is a cutout parallel to the z-axis through which matter as well as field 
can pass and leave the grid. Therefore, in these simulations no large-scale field can become 
attached to the black hole or fill the axial region as is seen in the GR simulations. Because 
the focus here is on the evolution of the disk away from the black hole and plunging region, 
we accept this limitation for the sake of concentrating computational power on the torus 
itself. 

The grid used in this study is designed to place as many zones as possible within the main 
body of the accretion disk. The inner boundary of the radial grid is set at Rin = 4:M. From 
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there it runs outward with constant AR to Ri, beyond which AR/R is set to a constant. 
This definition produces a logarithmically stretched grid throughout most of the domain 
while avoiding overly small AR values near the inner boundary. The z grid concentrates half 
of the total number of zones symmetrically around the equator using equal Az. This portion 
of the grid extends to z = ±zi from the equator. The remaining z zones are logarithmically 
stretched outward to boundaries that are well-removed from the initial torus. The grid 
uses evenly spaced zones that span one quarter of the full 2%. 

We choose a fiducial grid, Grid M ( "medium resolution" ) and vary the resolution relative 
to it. Grid M contains 256 radial zones with 48 zones inside of Ri = 20M. The outer radial 
boundary is at i? = 253M. The z grid uses 288 zones with 144 grid zones covering the range 
\z\ < Z\~ 5M. Using the density moment definition of the scale height, R' ~ 3.3M at the 
location of the pressure maximum, giving ~ 48 z zones per H' . Outside oi z — ±5M the 
grid is logarithmically stretched out to z — ±54M. The ^ grid uses 64 equally-spaced zones. 
In the inner disk (i.e., 6M < i? < 20M and \z\ < 5M), the cell aspect ratio AR/Az = 4.8 
and AR/{RA(f)) = 2.3-0.68 (from R^6M to R^ 20M). 

The simulations performed on Grid M will be contrasted with the results obtained using 
other grid resolutions. The lower resolution Grid L ( "low resolution" ) is half as well resolved 

in R and z. It consists of 128 radial zones with 24 equally spaced zones between R — AM 
and 20M, with AR increasing oc R from that radius out to R = 253M. At the pressure 
maximum, AR = 1.07M. The z grid has 144 zones, half of which are concentrated inside 
z = ±5M, where the minimum Az is 0.139M. The grid is the same as in Grid M. In this 
case, AR/Az = 4.8 as for Grid M, but AR/{RA(/)) = 4.5-1.4. 

Grid R ("high-R resolution") is the same as Grid M for z and </>, but the number of 
radial zones is increased to 816. Inside Ri = lOM, AR = 0.04M; outside that radius, 
AR/R = 0.004. The outer boundary of the radial grid is located at R = 142M. Its 
cell aspect ratios are AR/Az = 0.58-1.2 (once again, from R = 6M to R = 20M) and 
AR/{RA(f)) = 0.27-0.16. 

Two other grids are used to study the influence of the resolution. In Grid PL, the 
number of (p zones is reduced to 32. It uses 256 R zones, but redistributes them to decrease 
AR. For this grid, AR = 0.125M inside of = lOM, and AR/R = 0.0124M outside that 
point. Grid PH is the same as Grid M but with the number of zones increased to 128. 

We employ a variety of diagnostics to analyze the simulations. Certain quantities, such 
as total magnetic and kinetic energies, are computed and recorded every 10 timesteps. We 
also save binary data files and compute integrals of quantities over cylindrical shells at regular 
time intervals. 



- 9 - 



3. Shearing Boxes and Resolution Diagnostics 



To quantify the resolution effects in global simulations we begin with the shearing box 



( iHawley et al.lll995l ). a system where greater effective resolution can be employed and more 



extensive resolution studies are possible. Since shearing box models were first introduced, 



the initial " 


ield strength and geometry, and other thermodynamic factors (e.g.. Hawlev et al. 


1995. 


1996; 


Brandenbure et a 


1995; Stone et al. 


1996; Flemine et al. 2000: Miller & Stone 


2000; 


Sano & Inu 


;suka 


2001; ^ 


jano & Stone 2002; 


Sano et al. 2004; Fromane & Papaloizou 


2007; 


Simon et al. 


2009 


) The first shearing box simulations employed comparatively few total 



zones (e.g., 32 x 64 x 32 zones in x, y, z), but recently the number of zones used in simulations 
has substantially increased. 

Resolution studies using stratified shearing boxes that include the vertical component of 
gravity are the most relevant for comparison to global simulati ons. We have examined data 



from several recent stratified shearing box res olution studies (ISimon et al.ll201ll ; IShi et al. 



20 id ; iDavis et al.ll2010l ; iGuan fc Gammiell201l[ ) to see what, if any, general properties of the 
MRI turbulence might prove sensitive to resolution. Although shearing box simulations have 
found that the magnitude of a physical resistivity a nd viscosity as well as their ratios can have 
a sig nificant impact on the turbulence levels (e.g. Fromang et al. 2007; Ihesur fc Longaretti 
20071 ). at least for smaller magnetic Reynolds numbers ( lOishi fc MacLowll201ll ). nearly all 
global simulations done to date have used only gridscale dissipation. Therefore, we consider 
here only those shearing box simulations without explicit (physical) small-scale dissipation. 
The resolutions in these simulations range from 16 to 128 zones per scale height H. The 
Simon et al. and Davis et al. models have an isothermal equation of state and make us e 
of the Athena code with the HLLD-fiux solver (jStone et al.ll2008l ; IStone &: Gardinerll2010l ). 
but differ in regard to the initial magnetic field and box size. The Davis models use a 
zero- net vertical field initial condition and a box that is H : 4H : AH {x : y : z), while the 
Simon models have a net toroidal field overlaid with a poloidal field loop and a box that is 
2H : AH : 8H (x : y : z); in both cases, H = \/2cs/^. The Shi et al. models have the same 
initial field configuration as Simon, but use a Zeus code and include hea ting and radiative 
transport. Their box had side lengths 2H' : 8H' : 16H' (x:y:z). The iGuan &: Gammie 



(I2OIII ) models use a Zeus code, an isothermal equation of state, and focus on larger stratified 
boxes; their fiducial simulation h as a domain of 16H : 20H : lOH. T a ble U presents sorn e 
time -averaged data col l ected from ISimon et al.l (120111 ) , IShi et al.l (120101 ) , iDavis et al.l ( I2OIOI ) , 
and iGuan fc Gammid (120111 ) along with two unpublished 8 and 16 zones/if simulations 
(Simon, private communication). Note that the Shi et al. data use a time-averaged scale 
height because the temperature in their simulations varies as determined by thermal balance 
between dissipation of the turbulence and radiative cooling. 
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It is immediately clear that resolution can strongly affect the quantity most important 
to disk evolution, the magnetic stress. Results from four stratified shearing box simulations 
that use 64, 32, 16 and 8 zones per H for the first 150 orbits are shown in Figure [2J The 
plotted quantity is the volume-averaged stress parameter a, defined 

_ {{pvjvy- B^By/Ar:)) 
OL -— (^yj 



The 64 and 32 zone simulations are the 64Num and 32Num models of ISimon et al.l ( 120111 ): 
the 16 and 8 zone simulations are lower resolution versions of those simulations (Simon, 
private communication). All four simulations show initial growth to a peak that occurs at 
around 10 orbits in time. After this point, the 8 zone simulation (the lowest curve) clearly 
dies out. Large temporal variations characterize the other simulations. The stress in the 16 
zone run, for example, declines slowly for the first 50 orbits, but then rises by a factor of 3 
by orbit 60. The 32 and 64 zone runs maintain a time-averaged a that is consistent with the 
initial peak value. The 16 zone ru n also varies strong ly with time, but with a mean alpha 



that is less than in the 32 zone run. ISimon et al.l (120111 ) note that in unstratified simulations 



the time-averaged a changes by a larger amount in going from 16 to 32 zone resolution than 
in going from 32 to 64 zones per H (their Figure 1). In both the Davis et al. and Shi et al. 
resolution studies, the smallest number of cells per scale height was ^ 20-30, and there is 
little change in a when finer resolutions are employed, includir ig the best-resolved Day is et 



al. simulation, in which there are 128 cells per scale height. In iGuan fc Gammid ( 120111 ) the 
value of a went from 0.013 to 0.023 when the resolution of their fiducial model was doubled, 
a change from 13 to 26 zones per H in the vertical direction. 

The converged value of a in these particular stratified isothermal shearing boxes is of 
order ~ 0.02. There are, of course, additional effects beyond resolution that determine a. 
As discussed above, resistivity and viscosity can have a significant impact. At the same 
resolution, Simon et al. find consistently larger values of the stress compared to Davis et 
al., sometimes by close to a factor of t wo. The Simon et al. bo x is a factor of two larger 



in both the x and z dimensions. The iGuan &: Gammiel (l201l[ ) simulations also provide 



some evidence that a can be larger when larger domains are used. There is evidence from 
unstratified shearing boxes that taller boxes also promote stronger magnetic field (Stone, 
private communication). The Shi et al. simulations, in which the equation of state directly 
balances heating and radiative cooling, suggest that a may be somewhat larger when more 
realistic thermodynamics are employed. 

Even using the local shearing box, few simulations are carried out with resolutions as fine 
as 128 cells per scale height, so in practice one often asks the question, "In this particular 
simulation with only modest resolution, how close is the measured value of stress to the 
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numerically converged value?" Because high-resolution simulations can be very expensive in 
computer time, it is useful to define metrics of simulation quality that can be calibrated to 
a set of standard simulations and then applied to the data of a new simulation. In this way, 
how close that simulation comes to convergence may be estimated without the expense of 
additional, higher resolution simulations. In the remainder of this section, we discuss several 
such quantities. 



3.1. Convergence metric 7^1: 



The first such metric comes from the linear theory of the MRI. iNoble et al.l (l2010f ) used 
the vertical field characteristic MRI wavelength to compute a quality parameter defined 
by 

= Xmri/Az = (10) 

where Vaz is the z component of the Alfven speed. The characteristic wavelength Xmri 
is close to, but not precisely equal to, the fastest growing MRI wavelength. Wavelengths 
A < Xmri/\/S are stable, while all wavelengths A > Xmri are unstable, albeit w ith reduced 



growth rates oc (k • Vaz). On the basis of unstratified shearing box simulations, ISano et al. 



(120041 ) suggested that a value greater than 6 was required in order to achieve a linear 
growth rate close to the analytic prediction. Considering an isothermal thin disk with only 
vertical field in the initial condition, Xmri can be rewritten in terms of the plasma /3 by 
noting that (3 = 2pH'^VL^ / B"^ , and hence Xmri = 2nHP~^^'^. Thus, a value of Qz of ~ 10 
requires 1.6/3^^^ zones per H when the field is purely vertical; when the field has any other 
sort of geometry, /3 in this expression should be scaled by the fraction of the field energy in 
the vertical component, giving a zone total of 

iV. ^ 16 (/3/100)^/^ {{vl)/{v\^)Y^' (Qz/lO) (11) 

per scale height H. Because the fraction of the magnetic energy in vertical field is often only 
~ 0.01-0.1, the number of zones required for a given /3 increases by ~ 3-10. 

The second column of Tabled] shows the values of Qz, averaged over the midplane region 
{\z\ < 0.5H) for our stratified shearing box simulation sample. It is necessary to pick out 
the midplane region because \vaz\ generically increases sharply away from the midplane. 
Consequently, in these simulations in which the vertical resolution is uniform (unlike typical 
global simulations), Qz generally increases by 1-2 orders of magnitude from z = to 2; ~ 3H. 
These regions with better effective resolution can be important in maintaining the turbulence. 
By \z\ = 2-3H, /3 < 1, and the MRI is largely suppressed and the large values of Qz are 
less relevant. Comparing the Davis et al. series with the Simon series, we see that even 
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with 32 cells per scale height, the Sano et al. criterion is met only marginally. Although 
a increases dramatically when rises past a few, its dependence on resolution (in relative 
terms) appears to level out in the range 10 < < 20. 



3.2. Convergence metric ^^2: Qy 

Maintenance of poloidal field and turbulence requires non-axisymmetric motion. To 
estimate how well non-axisjmimetric stirring is described by the simulation, we can define a 
merit parameter Qy based on the toroidal field and the y grid zone size [Q^ and i?A0 for 
global simulations). The toroidal field MRI is nonaxisymmetric, and the linear properties of 
those nonaxisymmetric modes are somewhat different from those of the vertical field MRI. 
Although the nonaxisymmetric MRI modes depend on toroidal field, the presence of weak 
poloidal components can greatly increase the total a mplification of nonaxisy mmetric modes 



beyond what is predicted for a purely toroidal field (iBalbus h Hawleylll992l ). Like the case 
of vertical wavevectors, the maximum linear growth rate occurs for wavelengths comparable 
to the distance an Alfven wave travels in one orbit, but mode growth also depends on the 
radial wavelength, which evolves due to shear. Further, maximum gro wth also demands 



vertical wavenumbers much greater than H ^ (IBalbus fc Hawleylll992h . For shearing box 



simulations, the number of y zones required to achieve Qy ~ 10 is 

Ny ~ 64 (i//4) (/3/100)'/' (g,/10) (12) 

for a y direction spanning 4if. For toroidal modes in global simulations, Q^ = 2tcH/ (/3^/^i?A0), 
where P includes only the toroidal field component. To resolve linear growth of the toroidal 
MRI in a full 27r simulation requires 

- 1000 (0.1 R/H) (/3/100)^/2 (g^/io) (13) 

azimuthal cells. 

The simulations described in Table [H are nearly all well-resolved by the Qy criterion; 
this is one of the advantages of shearing boxes. Only in one case (the Simon 8 cells per 
H run) is Qy < 10. In both the Davis et al. and Simon et al. simulations, the cells are 
cubical. Because shear ensures that the azimuthal component of the magnetic field is much 
stronger than the vertical component, cell sizes that are too coarse to yield good vertical 
resolution can nonetheless be quite adequate to describe azimuthal behavior. However, in 
simulations whose grids are elongated in the azimuthal direction, Qy values will be smaller. 
In the better-resolved Davis and Simon simulations for example, Qy/Qz ~ 4 and one might 
expect that if /S.y//S.z ~ 4, Qy would be only comparable to Qz- 
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3.3. Convergence metric 7^3: a 



mag 



Whereas Qz and Qy derive from tlie critical linear wavelength of the MRI, the other 
diagnostics we have studied are more closely related to nonlinear development of the MHD 
turbulence. That is, they reflect how well a numerical calculation replicates macroscopic 
magnetic field properties related to the stress that are independent of discretization (i.e., 
convergence). 

The first of the nonlinear diagnostics is the ratio of the Maxwell stress to the magnetic 
pressure, defined amag = —'^BrB^/B'^. Although turbulent Reynolds stress also contributes 
to angular momentum transport at a level roughly 1/4 of the Maxwell stress, it is difficult 
to quantify in global simulations, so we do not include i t in our definition of amag- Even 



the earhest shearing box simulations (IHawley et al.lll995l ) found th at this quantity was re- 



markably constant from simulation to simulation. More recently, iBlackman et al.l (120081 ) 
examined a large sample of published unstratified shearing box results and found that, quite 
generally, a/3 cx amag is roughly constant, where a is the traditional constant of proportion- 
ality between stress and (gas) pressure; the combina tion al3 simply removes the gas pressure 



from consideration. Recomputing the results from iBlackman et al.l (120081 ) in terms of our 
definition of amag (i-e., without the Reynolds stress), we find that their derived values are 
more or less consistent with amag = 0.3-0.4 as seen in the shearing box simulations reported 
here, except for the one with only 8 cells per H. 

The values of amag presented in Table [1] are derived by taking the ratio of the Maxwell 
stress integrated over the central scale height to the similarly integrated magnetic pressure 
and then time-averaging; this is how amag has been determined in past shearing box sim- 
ulations. Other averaging procedures, such as averaging the local ratio rather than taking 
the ratio of the averages, give values that are generally somewhat smaller. The data show 
that at a gross level, amag and a are correlated: a very low value of amag corresponds to a 
very low value of a. In the lowest resolution simulation discussed here, with only 8 cells per 
vertical scale height, amag = 0.08. However, with an even modest improvement in resolution, 
both a and amag rise. 

The value of amag is relatively constant because in MRI turbulence and By are highly 
correlated. Both the background shear, which creates toroidal field out of radial, and the 
action of the MRI itself, which stretches out radial field as angular momentum is transferred 
between fluid elements, create this correlation. Time-averaging and volume-averaging over 
the central scale height, we find that the Pearson correlation coefficient C{Bx,By) is -0.73, 
-0.70 and -0.67 for the 32, 64, and 128 zones per H Davis runs, and -0.28, -0.75, -0.73, 
and -0.71 for the 8 through 64 zones per H Simon runs. In other words, once one is past 
a resolution threshold (between 8 and 16 zones per scale height), the correlation rapidly 
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achieves a value ^ —0.7. We have computed the correlation averaged over the [x, y) plane 
as a function of z. For \z\ < 2 it is consistently ~ —0.7, but approaches ze ro at higher 



altitudes. Two scale heights from the midplane is where the stress dies out (ISimon et al. 



20111 : iGuan fc Gammiell201l[ ). 



We can probe a bit deeper into the nature of amag by studying its probability distribution 
over the set of grid cells. Figure [3] shows a time-averaged distribution function for a^nag for 
the region within 2 scale heights of the equator for each of the four Simon runs. The low 
resolution distribution function peaks around 0; there is very little net stress remaining within 
the decaying turbulence, even while the azimuthal field persists. As resolution increases, the 
distribution shifts to higher values, as does the mean value. On a zone-to-zone basis, amag 
is correlated with Pmag in the sense that where the field is particularly strong, the ratio of 
stress to magnetic pressure is also particularly high. In order of increasing resolution, the 
mean amag values are 0.057, 0.275, 0.346 and 0.380. Thus, improving resolution also leads 
to greater correlation between and By, but the average a^ag saturates at ~ 0.4. 

Before leaving this topic, we note that the resolution-dependence of amag illustrates an 
important aspect of convergence-testing. If one compared its behavior in the 8 zone run 
with simulations having fewer cells per scale height, one might have concluded that amag is 
always ^ 1. In other words, low resolution simulations can be entirely blind to important 
effects, and convergence does not even begin until a resolution threshold is reached where 
that effect is at least minimally resolved. 



3.4. Convergence metrics #4 and #5: (Bl/By) and {Bl/Bl) 

The amag parameter depends, in part, on the relative magnitude of the poloidal and 
toroidal magnetic field components and in part on the degree of correlation between the 
radial and toroidal components. As we have already seen in our discussion of the Qz and Qy 
diagnostics, the fidelity with which a given simulation follows the linear growth of these two 
field components can be different. The same may be true of their nonlinear characteristics. 

To explore this question in a way that focuses on the separate components, independent 
of their correlation, we examined the time- and volume-averaged energy ratios (B^/By) 
and {BI/BD as functions of resolution. A relatively clear trend er nerges from the for mer. 



Figure m shows the time history of this ratio for the simulations of ISimon et al.l (120111 ). It 
increases from < 0.01 to 0.15-0.18 from the worst to the best resolved models. Despite the 
substantial time variability, at every stage the curves in this figure are well-separated from 
each other, demonstrating a systematic increase of this parameter with resolution. When 
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one looks at the full ensemble of shearing box simulations, the dependence on resolution 
becomes even more striking (Fig. |5]), with the value apparently leveling off near 0.2 at the 
highest resolution, when there are at least ~ 40 cells per H. On the other hand there seems 
to be no general trend for the ratio of vertical to radial field energies, {B'^/Bl); the values 
range from ~ 0.4-0.6 in these simulations. 

Moreover, as we will discuss in greater detail later, there is a strong correlation between 
the Q values and {Bl/By). The simulations achieving near-saturation values of {B^/By) 
(Simon64, Davis64, Davisl28, and ShiDBLE) all have Q, > 10 and Qy > 32. 



3.5. Summary: Stratified Shearing Box 



In principle, a resolution study is directly applicable only to a set of simulations using 
the same numerical scheme on the same problem. Only after cross-comparison of parallel 
resolution studies can different algorithms be calibrated relative to one another. Previous 
comparisons between the Athena and Zeus simulations have indicated that Athena has lower 
turbulence decay rates f or two-dime n sional shearing sheet sirn ulations compared to Zeus 
( IStone &: Gardinerll2010l : IStonell2009l ). IStone fc Gardinerl (120 lOl ) suggest that this is due to 
the use of third-order, rather than second-order spatial interpolation in Athena as well as 
the use of the HLLD-fiux solver. In the studies gathered he re, the contrast b e tween Zeus 
and Athena is less obvious: the Shi et al. simulations and the iGuan fc Gammid ( 1201 ll ) Zeus 
models have diagnostic values comparable to equivalently resolved Athena runs. 

Thus, from these simulations, we can conclude that for both Athena and Zeus simu- 
lations, stratified shearing boxes begin approaching convergence when resolution is around 
40 zones per H. The quantitative changes that result in increasing the number grid zones 
beyond this point are noticeable, but small, compared to the decrease in zone size, i.e. con- 
vergence is occuring at a rate faster than linear in Ax. Adequate resolution requires both 
Qz and Qy to be sufficiently large, but an especially large value of one can somewhat com- 
pensate for a smaller value of the other. When Qy > 20, Q^ > 10 suffices; however, when 
Qy is smaller, Qz > 15 is required. A ratio of radial to toroidal magnetic energy greater 
than > 0.15 and amag — 0.3-0.4 are signatures of well-developed MRI-driven magnetic tur- 
bulencej^ The ratio of the vertical to radial magnetic field energy, on the other hand, shows 
no particular trend with respect to resolution. 



^ These values are obtained from shearing boxes without a net vertical field or applied resistivity or viscos- 
ity, as appropriate for the global simulations to be discussed here. Further work is required to characterize 
shearing box turbulence in the presence of net vertical field and non-ideal MHD effects 
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It is likely that the importance of Qy stems from the essential role played by nonax- 
isymmetric processes in maintaining poloidal field energy and turbulence. Purely toroidal 
fields can support an active MRI-driven turbulence with relatively small vertical fields, but 
the vertical field MRI will necessarily generate toroidal field and nonaxisymmetric motions 
are essential to maintaining the poloidal field. 

Moving in the direction of reduced resolution, we note that even if the characteristic MRI 
wavelength is unresolved (small Q) longer wavelengths will be unstable, albeit with smaller 
growth rates. Thus, lower resolution simulations can still have MRI-induced turbulence and 
stress, but at correspondingly reduced levels. Such a state will be indicated by smaller relative 
values of the convergence metrics. Of course, the results from the 8 zone per H simulation 
show that there is a resolution limit beyond which turbulence cannot be sustained. 

4. Global Simulations 

In this section we describe the results from a set of pseudo- Newtonian global disk sim- 
ulations intended to investigate the influences of grid resolution and initial magnetic field 
strength and topology. What they show about progress toward numerical convergence will 
be discussed in the following section. 

Our parameter study is centered around a fiducial model, designated twoloop-1000- 
mr, a simulation with initial /3 = 1000, a two- loop initial magnetic field, and grid M, the 
medium resolution grid. Using the same initial condition, we have explored the effects of 
both increasing and decreasing the grid resolution. To study the effects of different initial 
magnetic field configurations, we have also used the M grid for a simulation with a single 
dipolar loop initial condition and the same initial /3, as well as both single and double dipolar 
loops with initial /3 = 100. The various models are listed in Table [2l 

4.1. The Fiducial Run 

The /3 value given for a torus is defined as the ratio of the total thermal to total magnetic 
pressure. Although the mean initial /3 = 1000, the /3 value at any given point within the 
torus varies. The minimum value is /3 = 283 and occurs along the radial field lines above 
and below the equator in the inner field loop; where the two loops meet at the pressure 
maximum, (5 = 412. For those regions of the torus where there is no magnetic field, /3 
is nominally infinite. The initial resolution parameter has three maxima, each located 
where the vertical field crosses the equator. From inside out, these maxima are Qz = 15, 29, 
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Table 1. Shearing Box Simulations 



Reference 


Zones/i? 


Qz 


Qy 


Bl/Bl 


Bl/Bl 


f- 


J-1 






a 


SimonS (unpub) 


8 


0.08 


1.7 


0.016 


1.0 





015 





08 





001 


Simonl6 (unpub) 


16 


2.0 


13. 


0.075 


0.53 





057 





30 





016 


Simon32 


32 


5.7 


27. 


0.13 


0.53 





072 





37 





025 


Simon64 


64 


11. 


44. 


0.17 


0.53 





056 





40 





020 


Davis32 


32 


4.5 


23. 


0.12 


0.41 





078 





33 





020 


Davis64 


64 


10. 


40. 


0.16 


0.47 





051 





36 





012 


Davisl28 


128 


26. 


98. 


0.18 


0.50 





053 





36 





018 


ShiSTD 


27 


4.8 


13. 


0.10 


0.65 





075 





27 





020 


ShiZ512 


53 


11. 


13. 


0.12 


0.56 





130 





30 





029 


ShiDBLE 


50 


15. 


32. 


0.15 


0.63 





098 





22 





029 


Guan stdl6 


12.8 


2.6 


15. 


0.07 


0.58 





035 





28 





013 


Guan sl6a 


25.6 


6.8 


34. 


0.12 


0.58 





057 





32 





023 



Table 2. Global Simulations 



Name 


Grid 


Type 


Initial /3 


Duration (M) 


twoloop-lOOO-mr 


M 


2- Loop 


1000 


66000 


twoloop- 100-mr 


M 


2- Loop 


100 


41000 


twoloop-lOOO-lr 


L 


2- Loop 


1000 


41000 


twoloop-lOOO-hr 


R 


2-Loop 


1000 


34000 


twoloop-lOOO-mlp 


PL 


2-Loop 


1000 


40000 


twoloop-lOOO-mhp 


PH 


2-Loop 


1000 


25000 


oneloop-lOO-Ir 


L 


Dipole 1-Loop 


100 


40000 


oncloop-lOO-mr 


M 


Dipolc 1-Loop 


100 


39000 


oneloop- 1000- mr 


M 


Dipole 1-Loop 


1000 


38000 
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and 20; the linear growth phase of the MRI should be well-resolved. 

The fiducial model was run for 6.6 x lO^M in time, corresponding to 54 orbits at the 
initial pressure maximum, and over 1000 orbits at the location of the ISCO, R = 6M. This 
evolution time is longer than previous global simulations, and permits an examination of 
the long term behavior of the accretion fiow and MRI-driven turbulence at these resolutions. 
The evolution proceeds in several stages as identified by the evolution of the total magnetic 
energy (Fig. E]). First there is a period of exponential growth in the magnetic field. This is 
rather brief; hj t = lOOOM (about 1.3 orbits ai R = 24M), significant radial field has already 
grown to supplement the vertical field in the innermost part of the inner loop. This behavior 
is characteristic of the vertical field MRI. The total integrated radial field energy has doubled 
by this time; the vertical field energy doubles by t = 1500M. At first, the toroidal field grows 
as the MRI creates radial field and the background shear stretches radial field into toroidal; 
later azimuthal MRI modes also contribute to toroidal field amplification. By t = 2000M the 
total toroidal field energy is 10 times as large as the total poloidal field energy, and significant 
accretion (mass fiow off the inner radial boundary) has begun. Global field energy continues 
to rise until it peaks at t ~ 8500M. It then declines until t = 1.5 x lO^M, after which the 
volume-integrated energies of the different magnetic field components vary slightly around 
a slowly declining trend line. We will refer to the period between 2000M and 1.5 x lO^M as 
the "initial peak" ; the remainder of the simulation we call the "quasi-steady state" period. 



4.1.1. Inflow equilibrium and the quasi- steady state 

Nearly every global disk simulation begins with a finite mass on the grid. Because 
accretion of some mass entails transfer of its angular momentum to other mass, part of the 
matter in the simulation must move outward as other mass moves inward. Consequently, at 
most only a portion of the disk can actually be in a state of infiow. At best, therefore, infiow 
equilibrium can be established only within some radius. 

Having identified the largest radius within which infiow equilibrium can be sought, one 
might define infiow equilibrium to be a state in which the mass accretion rate is constant 
as a function of radius. The problem with this definition is that accretion is driven by the 
fundamentally chaotic process of MHD turbulence. The accretion rate at a specific radius 
must therefore always be highly variable in time. One solution is time-averaging; regions 
of infiow equilibrium would then be those ranges of radius over which the time-averaged 
accretion rate is constant. 
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A closely related procedure is to make use of the equation of mass conservation 

Here S is the surface mass density and M the accretion rate integrated over the cylindrical 
surface at radius R. Clearly, wherever M is independent of R, E is constant in time. 
Thus, one could also test S(i?, t) for time-steadiness across the radial range of interest. 
Equivalently, one could check that M(< R,t) = J dR' 2'kR'T.{R' ,t) is time-steady. This 
alternative has the conceptual advantage that it focuses squarely on the primary matter of 
interest: that the mass distribution in the disk does not change secularly over time. It also 
has the technical advantage that t) changes more slowly than M{R,t), so results taken 
from simulation data are less subject to noise fluctuations. 

Figure [7] shows M(< R,t) for the fiducial run. The disk mass rises quickly during the 
first lO^M in time and then levels off. After that time, there is a slow secular diminution 
in the mass of the disk within R = 20M, but its characteristic timescale is quite long, 
~ 5 X lO^M at R= lOM, ~ at i?, = 20M. Because the time to drain and replenish the 
region inside R = 20M is M/M ~ 1.5 x lO^M, this secular trend is quite slow compared 
to the characteristic mass equilibration time. Despite the overall steadiness of the radial 
mass profile, there are also shorter timescale fluctuations that become progressively larger 
in fractional terms at smaller radii, reaching ~ 50% at i? = lOM. These reflect, of course, 
the continuing large amplitude fluctuations in the mass accretion rate both as a function of 
time and of radius. The slow diminution in inner disk mass reflects the declining trend in 
the magnetic energy, which leads to a parallel fall in the mass accretion rate. 

Another way to test for inflow equilibrium is to compare the simulation time with an 
estimate of the characteristic inflow time. We can compute an average inflow velocity from 
simulation data by taking 

^"^(''^^ = jpRd4>dz ■ ^''^ 

These values are averaged over time to remove the ever-present fluctuations. The accretion 
time from radius R is then tm(-R) = dR! / {vr{R')) . Computing this for the fiducial run 
at i? = 20 we obtain 1.2 x lO'^M, consistent with the estimate above based on the average 
accretion rate and the disk mass interior io R = 20M. 

How does this result compare to an estimate obtained from steady-state disk theory? 
In the steady-state limit the equation of angular momentum conservation is 

Wr, = ^{1-3.I3) (16) 
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where Wr^ is the vertically averaged R(f) component of the stress tensor, and j and j* are 
the specific angular momentum at R and the angular momentum accreted per unit mass. 
The j^, term determines the net fiux of specific angular momentum; traditionally, it has been 
set to the angular momentura at th e ISCO on the assumption that stresses cease there. 
Following IShakura fc SunyaevI (119731 ) . we write the vertically- integrated stress in units of 
the vertically-integrated pressure and assume there is a single temperature throughout the 
fiow; i.e., we set Wr^I) = aSc^. With those assumptions, we obtain for the steady-state mean 
infall velocity 

Agreement between the infiow velocities from ( fT5|) and f|T7|) would indicate that the ob- 
served mass accretion is in approximate steady state at a rate consistent with the angular 
momentum transport produced by the observed stress. 

We test this for the fiducial run for a time average over 10*^M beginning at t = 4 x lO^M. 
We compute a vertically and time-averaged Maxwell stress, pressure, and density to obtain a 
and c^. The average value of a thus obtained for the Maxwell stress is 0.017. To account for 
the Reynolds stress (not directly measured in global simulations), we increase the value of a 
in f ll7p by 25%, a value consistent with results from shearing box simulations. The velocity 
derived from the steady state disk model (with = jisco) and the velocity obtained directly 
from the accretion rate and the mass are compared in Figure El The agreement is good out 
to just beyond 30M. This match is consistent with the range over which the time-averaged 
M is constant with radius, the computed at 30M which is 4.5 x lO'^M, and our estimates 
based on the evolution of M(< i?, t) as described above. The two curves deviate outside 
of i? ~ 32M because there the infall time begins to exceed the simulation time. At small 
radius the curves deviate because the stress does not go to zero at the ISCO; reducing to 
0.985^75(70 brings the curves into line right down to the ISCO. 



4.1.2. Comparison with shearing box results 

The initial transient and quasi-steady periods are seen in both global and shearing box 
simulations. Figure [9] plots the total magnetic field energy, normalized to its peak value, as 
a function of time in units of orbits at the initial pressure maximum. Overlaid on this is the 
evolution of the total magnetic field energy in the stratified shearing box simulation with 16 
zones per H, again normalized to the initial peak value. Because the early evolution of the 
global simulation is relatively local, dominated by MRI growth in the confined region where 
its growth rate is greatest (i.e., the inner rings of the initial torus), it is perhaps not surprising 
that the initial evolutions are similar. More interesting, perhaps, is the similarity between 
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the two models during the subsequent quasi-steady state phase. Between 20-50 orbits in 
each run, the ratios of the total average radial to toroidal magnetic field energy, {Bj^/B'^), 
are 0.068 and 0.070 for the shearing box and the fiducial model, respectively. This similarity 
suggests that the effective resolution in the fiducial global simulation, at 14-60 zones per 
H depending on radial location, is comparable to the 16 zones per H used in the shearing 
box. A comparison with Fig. [2] shows that higher resolution shearing box simulations have 
a smaller decline in stress immediately after the initial peak. The 16 zone shearing box 
simulation, however, sees a regrowth of field energy beyond orbit 50 that is not seen in the 
global simulation. The reason for this is uncertain. It could be due to the better azimuthal 
resolution in the shearing box, possibly by better capturing a (nonaxisymmetric) dynamo 
process. Another possibility is that the observed field regrowth is due to other properties 
associated with the shearing box, e.g., restricted box size, shearing-periodic boundaries, and 
overall symmetry. 

Among the differences between shearing box and global simulations is the latter's large 
dynamic range in Q, which is important because the MRI e-folds at a rate ~ As a result, a 
single spatially-averaged value is not very useful for comparison with shearing boxes. Instead, 
we show figures of the principal diagnostic quantities as a function of radius, averaged in 
time and azimuthally- and vertically- averaged weighted by density (Figs. [TO] and [TTl) . In each 
case, the averaging interval was chosen to be the longest time-span covering the quasi-steady 
epoch for all simulations shown in a given figure; consequently, the averaging periods for 
the different figures are different. Likewise in each case, the radial range was restricted to 
the region defined as the "inner disk," QM < R < 20M in order to focus attention on the 
portion of the simulation most closely resembling a statistically time-steady accretion fiow. 

Data from the fiducial run appear in both of these figures. Its density- weighted (Qz) 
drops steadily inward, from ~ 8 at i? = 20M to < 2 at i? = 6M. By contrast, its (Q^jj) 
value varies only slightly with radius, remaining near ~ 8-10 throughout the inner disk. 
Both behaviors — the strong dependence of (Qz) on radius and the near-constancy in radius 
of (Qfj)) — are characteristic of all the simulations. {B\/ B'^ displays a radial profile very 
similar to that of {Qz)i falling from ~ 0.08 in the outer part of the inner disk to < 0.04 just 
outside the ISCO at R = 6M. 

Thus, even with a poloidal cell-count of 256 x 288, the fiducial run is, at best, marginally 
resolved according to both the Qz and criteria. There is adequate resolution to describe 
linear MRI growth of poloidal perturbations only near R = 20M; nowhere in the inner disk 
is it well enough resolved to describe poloidal nonlinear behavior properly. The azimuthal 
resolution is no better: (Q^) never reaches the ~ 20 level indicated by shearing box simu- 
lations. Similarly, the value of the nonlinear criterion {Bj^/ BV) is at most only about half 
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what the shearing box simulations suggest is characteristic of well-resolved turbulence. It is 
near to the value seen in the marginally resolved 16 zones per H shearing box. 

The radial gradient in (Qz) is an illustration of how grid-design can interact with physics. 
The definition of this quantity is 2tt\vaz\/ {^^z). In the cylindrical coordinates used here, 
is independent of radius. Consequently, Qz oc vaz{R)R^^'^] unless vaz rises rapidly inward, 
poloidal resolution quality falls toward small radius. Here the average vaz goes roughly like 
inner disk so {Qz) oc R. By contrast, oc va4,R^^'^ / ^(t>i leading to its much 
weaker dependence on R. The end- result of these different dependencies on radius is that, 
unless Qz is very large at larger radii, the linear growth rate of axisymmetric MRI modes 
will be reduced in the inner disk even while non-axisymmetric modes continue to grow at 
their correct (albeit slower) rate. 

The radial gradient in (Qz) also provides a finely-sampled measure of convergence prop- 
erties. At each radial cell, these global simulations sample a different effective resolution. 
As we will discuss in § [5l this effect allows us a much more quantitative measurement of the 
convergence rate than we would otherwise be able to obtain. 

4.2. Different grid resolutions 

Table [3] lists time- and volume-averaged values for various parameters, both diagnostic 
and physical, as well as the time-averaged accretion rate through the inner R boundary 
as a fraction of the initial torus mass. The time averages are taken over the steady state 
period for each simulation; the volume averages are limited to the inner disk body. For 
those parameters with potentially significant radial gradients, a range of values is given: the 
first number refers to the value at i? = 6M, the second to R = 20M. The scale height 
H is defined by a time- and density- weighted mean of ^/2cs/^^■ In these simulations, in 
which total energy is not conserved, Cg scales with radius roughly oc R~^^'^, so that H is 
approximately oc R. In most cases, amag increases outward, but simulation twoloop-lOOO-hr 
is an exception: in that case, amag decreases slightly toward larger radius. 

4-2.1. Azimuthal Resolution 

We next consider how the results seen with the fiducial model change as the resolution 
is altered. We begin with a study of the resolution (Fig. [TTj) . Model twoloop-lOOO-mlp 
uses grid PL, which decreases the number of zones to 32 while gaining a modest increase in 
radial resolution due to bringing in the outer boundary somewhat. Twoloop-lOOO-mhp uses 
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grid PH, which is the same as the M grid, but with an increase in the number of zones to 
128. As the orbital speed generally sets the Courant limit in these simulations, this high-0 
resolution model is more costly to evolve. 

Figure [12] shows the evolution of the total poloidal magnetic energy as a function of 
time for the three runs. For the 32 zone run, exponential growth ends earlier compared to 
the other simulations, with a more gradual climb to a peak value. The 128 zone run leaves 
the rapid growth stage earlier than the 64 zone run, but when it levels off after declining 
from its magnetic energy peak, it does so at a higher level. Thus, the sustained field energies 
after the initial peak are clearly separated by resolution. 

As shown in Figure [TTb . for R > 12-15M, (Qcf,) increases by roughly a factor of 2 
for each factor of 2 improvement in resolution. Notably, however, both poloidal indicators, 
(Qz) and {Bj^/B"^), also respond positively to improvement in toroidal resolution at larger 
radii within the inner disk. Like (Q^), the response to improved resolution is even stronger 
between the M and PH grids than between PL and M. In other words, improved azimuthal 
resolution leads to stronger poloidal field, even in the absence of improved poloidal resolution, 
and this effect strengthens with finer azimuthal resolution. 

4-2.2. Poloidal Resolution 

Figures [13] and [10] show the diagnostics for one-loop and two-loop initial field config- 
urations with a range of resolutions. In the one-loop factor of 2 improvement in 
poloidal resolution yields an equal factor in {Bj^/ Bp and a factor of 3 improvement in (Qz). 
Moreover, even with a grid of 256 x 64 x 288 zones, these simulations remain unresolved. The 
best (Qz) was only ~ 10, while it fell below 6 inside R ~ 10; {Bj^/B'^ hovered around 0.1 
outside R ~ 15M, falling to 0.06-0.07 near the ISCO, whereas the shearing box simulations 
indicate a converged value of almost 0.2. 

Within the two-loop series of simulations, we find that the move from the L poloidal 
grid to the M yielded the same level of improvement as in the one-loop factor of 2 

in {Bj^/B'p and a factor of 3 in (Qz). However, progress stalls when only the radial zone 
size is reduced, as in the R grid. Outside R ~ 10-12, the finer radial resolution actually led 
to a deterioration in both {Bj^/ B'p and (Qz)- We speculate that this may be due to the 
highly-elongated cells created in this grid by its finer radial resolution. In contrast to the M 
grid, in which RA(j)/AR ~ 0.4-1.4, that ratio is 4-6 in the R grid, increasing toward larger 
radii. 

Although not illustrated in the figures, much the same can be said about the effect of 



Table 3. Resolution Comparisons 



Name 


R,4 


), z zones 


M 






Qz 


Q 









Cells per H 


twoloop-lOOO-mr 


256 > 


< 64 X 288 


2.71 X 


10" 


6 


0.24-0.33 


1.6-7.5 


7.9- 


-10. 


0.040-0.077 


0.0088-0.018 


11-60 


twoloop- 1 00- mr 


256 > 


< 64 X 288 


2.51 X 


10- 


6 


0.29-0.37 


3.0-9.0 


12.- 


-11. 


0.054-0.092 


0.012-0.023 


14-65 


twoloop-lOOO-mlp 


256 > 


< 32 X 288 


1.99 X 


10" 


6 


0.25-0.24 


4.2-4.5 


8.5- 


-6.0 


0.059-0.046 


0.016-0.012 


8-49 


twoloop- 1000- mhp 


256 X 128 X 288 


4.40 X 


10" 


6 


0.23-0.41 


2.1-11. 


22.- 


-22. 


0.045-0.13 


0.0079-0.035 


13-57 


twoloop- 1000-lr 


128 X 64 X 144 


0.94 X 


10" 


6 


0.095-0.25 


0.46-2.3 


5.7-9.4 


0.019-0.050 


0.0055-0.0081 


8-30 


twoloop- 1000-hr 


816 X 64 X 288 


1.89 X 


10- 


6 


0.31-0.29 


3.8-4.9 


13.- 


-11. 


0.064-0.057 


0.021-0.022 


10-46 


oneloop- 1000-mr 


256 X 64 X 288 


3.85 X 


10" 


-6 


0.29-0.35 


2.8-9.1 


12.- 


-11. 


0.059-0.086 


0.012-0.024 


14-60 


oneloop-lOO-lr 


128 > 


< 64 X 144 


1.43 X 


10- 


6 


0.13-0.29 


0.42-3.3 


5.6- 


-12. 


0.017-0.065 


0.0053-0.013 


9-36 


oneloop-lOO-mr 


256 X 64 X 288 


7.92 X 


10- 


6 


0.30-0.38 


3.0-10. 


13.- 


-12. 


0.063-0.10 


0.010-0.033 


13-76 
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resolution on amag- When the magnetic geometry is held fixed and only resolution changes, 
finer resolution in both the poloidal and azimuthal grid lead to larger values of amag- 

Overall, it appears that improving radial resolution without also improving either ver- 
tical or (perhaps especially) azimuthal resolution does not yield significant progress toward 
convergence. As we have seen, the linear indicators {Qz and Q^) do not grow (if anything, 
they fell in this case), nor do nonlinear indicators such as (Bj^/B'i) improve. 



4.3. Different initial magnetic fields 

All simulations must begin from some particular magnetic field. The problem, of course, 
is that we can only guess at what might actually occur in Nature. We have no way to say 
what geometric structure it might have, yet numerous variations are imaginable. The field 
may or may not have net fiux; even if it has no net fiux, it may be predominantly toroidal 
or poloidal; and numerous sorts of poloidal configurations are possible. Even granted the 
field geometry, we must still specify its initial strength; so long as it corresponds to a plasma 
/3 ^ 1, we have no guidance in this respect either. The most we can hope for is that in the 
end the resulting steady state accretion fiow is independent of our arbitrary choices. 

We begin with the question of sensitivity to initial field strength with a comparison of 
oneloop-lOO-mr and oneloop-lOOO-mr. Stronger fields, of course, will have an initial advan- 
tage of larger Q values at a given resolution. In the quasi-steady state period, as shown 
by Figure [131 ^ factor of 10 in initial magnetic field pressure makes a consistent, but very 
small difference: the stronger initial field in general has larger values of the diagnostics, but 
only by ~ 10%. There is a greater contrast between twoloop-lOO-mr and twoloop-lOOO-mr, 
generally ~ 20-30% and in some places larger. Like the one-loop case, the simulation with 
the initially stronger field retains that advantage. 

The main effect of a stronger initial field is in the initial evolution; it creates stronger 
accretion early on and fills the inner disk with matter at a faster rate. This happens for two 
reasons, neither of which has much long-term effect. The first is that the stresses driving the 
initial infiow are not created by correlated MHD turbulence, but by the shear that transforms 
the initial poloidal field into toroidal. The second is that in the initial state, (Qz) is smaller 
if the field is weaker, depressing the linear MRl growth rate in some locations. In the end 
how the mass distribution reaches its steady-state is not important, so long as one evolves 
long enough to reach the steady-state. 



We next turn to the question of field geometry. iBeckwith et al.l (l2008al ) demonstrated 
that the magnitude of the axial funnel field attached to the black hole, capable of powering 
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a jet when the black hole spins, depends strongly on whether the vertical component of the 
initial magnetic field in the simulation has a consistent sense over a wide span of radii — 
the jet luminosity is substantial only when it does. They also showed that the accretion 
flows in the cases they considered depended much less on the field structure, although the 
stress levels and accretion rates are noticeably lower when the initial field is purely toroidal. 
Here we take a closer look at the accretion fiows associated with two specific initial field 
configurations: one initial dipole loop versus two. 

The first point is that the two-loop configuration does not persist throughout the simu- 
lation. Distinguishing the loops by their sign of A^, the azimuthal component of the vector 
potential, we find that in the fiducial model the disk inside R = 20M is effectively dom- 
inated by the inner loop until t ~ lO^M. By the end of that period, much of the inner 
loop's fiux has been accreted through the inner boundary, and the remainder has risen to 
high altitude and left the disk proper. For the next ~ lO^M, the main body of the disk is 
entirely dominated by the outer loop. Between ~ 2.2 x lO'^M and ~ 3.2 x lO^M, enough 
inner loop fiux has settled into the disk that it is very mixed. During this period, there is 
significant reconnection between oppositely directed field loops. For the remainder of the 
simulation (until t = 6.6 x lO^M), the disk is once again dominated by poloidal field loops 
having the sense of the outer of the two initial loops. In other words, with the exception of 
an episode that lasts only ~ 1/6 of the simulation, there is little to distinguish the poloidal 
field geometry in the inner disk of a two-loop simulation from one that begin with only a 
single set of poloidal loops. In both cases, there is a consistent sign of A^j, (i.e., consistent 
sense of field circulation) in the inner disk, but with a highly complex structure (the more 
fully-developed the turbulence, the more complex, of course). 

Nonetheless, our time-averaged diagnostics do show interesting contrasts between the 
one-loop and two-loop results (Figs. [T3] and [10]). At low-resolution, the (Q^) diagnostic 
improves by ~ 50% from two-loop to one-loop at all radii. At medium-resolution, (Q^) in 
the one- loop case is only about 10% greater in the disk body (e.g., R = 20M) than for the 
two-loop case, but that advantage widens toward smaller radii, rising to almost a factor of 
2 by i? = lOM. Very similar contrasts are found for {Bj^/B'^). Although the magnitude of 
the change is smaller, the trend is the same for amag- 

One might reasonably ask why these contrasts occur, given our assertion that the inner 
disk magnetic field at any single time generally has a uniform sign of A^. A clue comes from 
looking at their time-variation. During the period 1-2 x lO'^M, the diagnostic quantities 
in twoloop-lOOO-mr are, in fact, very similar to those of oneloop-lOOO-mr. Their decline at 
small radii takes place only after 2 x lO'^M, the time when oppositely-directed field loops 
reenter the disk. That re-entry is extremely irregular and creates gradients in on very 
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short lengthscales. For the same gridscale, the larger amount of short lengthscale fluctuation 
power leads to a faster numerical reconnection rate. This, in turn, weakens the magnetic 
field to the point that the MRI is unable to rebuild the field strength, even after the period 
of enhanced reconnection is over. Throughout this later period Qz is only ~ 2-8, insufficient 
to drive the MRI at its full rate. 

Thus, the two-loop initial condition creates — through a rather recondite mechanism — 
considerable short length-scale power that requires a finer resolution grid to treat properly. 
In principle, a weakened magnetic field (and therefore smaller values of {Qz) and {Q^) for a 
fixed gridscale) might be a physical consequence. However, there is clear evidence that this 
really is a resolution effect. As we have just pointed out, the prevailing Qz after the strong 
reconnection epoch is inadequate to drive linear MRI growth, much less support nonlinear 
poloidal dynamics. The low value of {B\/ B"^) further supports our conclusion that the 
turbulence is less than fully-developed at late times in the two-loop simulations. 



5. Interpretation 

The principal goal of this paper is to develop from highly resolved shearing box models 
a set of diagnostics that can be applied to gauge how well global simulations represent the 
quasi- stationary behavior of a disk obeying the same physics as in the simulation. In this 
section we will discuss how to use the diagnostics we have identified, making use of their 
application to the simulations described in the previous section. 

One might have supposed that the most natural quantity to study as a gauge of progress 
toward convergence is the stress, the physical quantity of central interest to the entire subject. 
However, there are several good reasons not to use it. The first is that it is not dimensionless. 
A similar statement can be made about the accretion rate. It can be useful to compare 
quantities such as stress and accretion rate between a set of simulations with the same 
initial conditions, but a comparison of any dimensional value across different simulations 
would need calibration in terms of other quantities. 

The stress is often quoted in units of the gas pressure, giving it an apparently "natural" 
dimensionless form, but it is not at all clear that this form is appropriate. We do not know 
whether the pressure is the only variable influencing the str ess, or whether the dependence 



is linear. Fluctuation timing correlations (IHirose et al.ll2009l ) suggest that the relationship is 



more complicated than commonly thought. In addition, even if stress were dependent only 
on the pressure and exactly proportional to it, the great majority of accretion simulations 
treat thermodynamics with very crude approximations, making the value of the pressure 



- 28 - 



both questionable and difficult to compare from one simulation to the next. For all these 
reasons, therefore, we prefer to use diagnostics that are dimensionless ratios of quantities 
referring only to magnetic properties. 

5.1. Relationship between the different diagnostics 

Each of our diagnostics has a slightly different standing. (Qz) and (Qc/.), for example, 
depend directly on cell sizes in the grid and can be meaningfully evaluated on the initial data 
whenever it includes either vertical magnetic field (for {Qz)) or toroidal field (for {Q4,)). In 
that sense, they are predictive measures of the abihty of the simulation to describe accurately 
the linear growth of, respectively, the axisymmetric and non-axisymmetric branches of the 
MRI. At the same time, however, both are also simulation products because their values 
change as the field strength evolves; in that sense, they also serve as diagnostics for the fully 
developed MRI-driven turbulence. 

The quantity {Bj^/B'^) has little significance in the initial state; it is simply a function 
of the initial magnetic field. It is not even defined if the initial field is entirely poloidal. 
Thus, it is purely a product of the simulation; the same is true of amag- 

The definitions of the Q parameters have an explicit dependence on grid size; their 
values can never converge with increased resolution. In contrast, (Bj^/B^) and amag measure 
physical quantities and do not have an explicit dependence on resolution in their definitions. 
In principle, they could be measured in a laboratory experiment. Thus, they signal progress 
of the simulation toward achieving the true physical state of the MHD turbulence, no matter 
what resolution or algorithm has been used. 

Despite these contrasts in status, there is a tight relationship between these different 
diagnostics of resolution quality. This relationship is illustrated in Figure [TH Although the 
filled circles in this figure are drawn from 11 different times at all radial cells in the inner 
disk in each of 3 different simulations (twoloop-lOOO-lr, twoloop-lOOO-mr, and twoloop-1000- 
mhp) that differ in both poloidal and azimuthal griding, they closely follow a single track: 
(ByBl) ~ 0.01 ((5^)°-^^ for (Qz) < 10-12, flattening out at larger (Q^), where it rises slowly 
to almost 0.2. The crosses, squares, and diamond represent the time-averaged values in the 
midplane region of the different stratified shearing box simulations; the different symbols 
distinguish different ranges of (Qy) (the crosses have similar {Qy) to (Q^) in the global 
simulations; the others have larger (Qy))- The shearing box simulations with (Qy) similar 
to the (Q^) values of the global simulations lie very close to the track defined by the global 
simulations. Consonant with our finding that better azimuthal resolution also strengthens 
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poloidal field, the shearing box results with larger (Qy) lie somewhat above the track; that is, 
larger {Qy) can compensate to a certain degree for smaller (Qz)- Taken as a whole, this figure 
confirms that for global simulations as well as for shearing boxes, convergence in (Bj^/B'p 
does not begin until {Q^) > 10-15 for typical global simulation azimuthal resolution (i.e., 
(Q^) ~ 10); larger (Qy) (> 25) can relax this requirement slightly. 



5.2. Importance of ?//0-resoIution 



In most of the MRI literature hitherto, discussion of resolution cri teria has generall y 
focused on poloidal griding, particularly the criterion established by ISano et al.l (12004 ). 
Less attention has been paid to the resolution. It is generally understood that in multi- 
dimensional simulations one should avoid cell aspect ratios that are too far from unity. 
Accuracy in directionally split algorithms requires that the truncation errors in the different 
spatial directions not be too different. In any case, if the truncation error in one spatial 
direction remains large due to a lack of resolution, improvements in resolution in other 
dimensions will be of limited value due to the dominance of that error. 

In the context of disk dynamics, including shearing box simulations, experience suggests 
that even though it is best for AR/Az ~ 1, ratios RA(f)/AR as large as ~ 4 are acceptable 
because orbital shear tends to stretch out features in the azimuthal direction. Our findings 
here suggest that larger ratios, such as are found in twoloop-betalOOO-hr, tend to weaken the 
development of MHD turbulence even though the unfavorable ratio was created by improved 
radial resolution. 

To our knowledge, the comparisons presented here are the first to underline the impor- 
tance of adequate (Q,^) . Our analysis indicates that sufficiently large (Qip) is necessary for the 
proper development of poloidal properties such as {B\/ B"^) and maintenance of a large {Qz)- 
This is, perhaps, not so surprising as it is well-understood that non-axisymmetric motions 
are essential for maintaining poloidal field in the face of dissipation (the "anti-dynamo" 
theorem). Ev en the earliest stratified shearing box simulations found evidence for a dy- 
namo process (IBrandenburg et al.lll995l ). and recent simulations have shown explicitly that 
magnetic field evolution can be empiricall y modeled with a simple a-Vt dynamo equation 
(IGuan fc Gammidl201ll : ISimon et al.ll201ll ). Although we have yet to develop a detailed un- 
derstanding of the mechanism by which the MRI creates the a-effect that generates poloidal 
field from toroidal, it clearly requires adequate resolution to be effective. Future simulations 
should therefore be graded on their values of {Q^) as well as {Qz)- 



Finally, one should bear in mind that resolution might well be much more important 
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in global simulations compared to local. The shearing box is in a local co-rotating frame 
whereas the global simulation must deal with the full orbital velocity. If one assumes that the 
numerical dissipation level in a code is proportional to a diffusion coefficient, e.g., Ax^/At, 
then the diffusion will be roughly proportional to the fastest velocity times Ax because 
At ~ f^/Ax, where is that fastest velocity. For a global simulation, the orbital velocity 
is Rfl oc {R/H)cs. Because for a shearing box the fastest velocity is < Cs, one might well 
require R/H more grid zones in the direction compared to a shearing box in order to 
have similarly low levels of diffusion. For global thin disk simulations this requirement is 
particularly problematic. 



5.3. Connection with stress 

It would seem natural that stronger magnetic fields would lead to larger stress. In fact, 
the correspondence between the hierarchy we see in our diagnostics {{Qz), {Q<j,), 
c^mag) IS rcproduccd closely in stress levels. Figure [T^ plots the stress as a function of radius 
and can be compared with the previous figures showing diagnostic quantities versus radius. 
As the diagnostics indicate better resolution, the stress rises. Moreover, when poloidal 
resolution is improved and the general level of stress rises, so does the ratio between stress 
near the ISCO and stress in the disk body. Although not illustrated in the figure, this trend 
also applies to the radial stress profile at particular times during an individual simulation. 
Higher stress levels at larger radius generally correspond to higher stress throughout the 
flow, even if the stress and the Q values decrease moving in. The ISCO-region stress appears 
even more sensitive to the prevailing stress than the stress in other regions of the inner disk. 

Because even relatively low resolution simulations nevertheless show some stress due 
to MRI-turbulence, the mere presence of stress and accretion are not themselves indicative 
of convergence in any sense. Weak flelds remain unstable at longer (better resolved) wave- 
lengths, but at greatly reduced growth rates. As we have tried to make clear, not even the 
best-resolved global simulations shown here appear to be near the diagnostic levels identifled 
as adequate in the shearing box models. Even if some regions satisfy the resolution criteria, 
their variation with radius means they are not satisfled globally. We expect that still better- 
resolved simulations can be expected to show both higher stress levels overall, and a greater 
ratio of ISCO-region stress to stress in the disk body. 
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5.4. Effects of the initial magnetic configuration 

Because the MRI grows exponentiaUy when the initial field is weak (/3 ^ 1), it has long 
been hoped that the saturated state of this instability would be independent of the initial 
field intensity. The stre ss in unstratified shearing box simula tions depends on the net value 



of the initial field (e.g., iHawley et al.lll995l : ISano et al.ll2004l ). but this probably refiects the 



constraints imposed by the shearing box boundary conditions, which preserve net magnetic 
fiux. In less-constrained global simulations one might expect the field to evolve to some 
self-consistent average magnetic energy with /3 somewhere below equipartition, independent 
of the field's initial st rength, even though local regions may have persistent net vertical fiux 



(jSorathia et al.ll2010l ). Although we have not made extensive comparisons of simulations 
with different initial values of /3, those we have examined are consistent with this hope. In 
the quasi-steady state period there are only modest differences between our 1- and 2-loop 
simulations with initial /3 = 100 and /3 = 1000, for example. 

Field geometry can matter, however, in the way in which it interacts with grid resolution. 
Intuitively, one might suppose that a correct discretized description of any field requires a 
grid with a characteristic scale significantly smaller than the field's characteristic lengthscale 
of variation. Applying that argument in the present context suggests that initial fields with 
more small-scale structure require finer grids to describe accurately. The contrast between 
our results for the 1-loop and 2-loop simulations supports this suggestion, although, as we 
have explained, the specific way in which the small-scale structure is created is rather more 
complicated than one might have guessed. Introducing oppositely-directed field-loops creates 
more opportunities for reconnection, and this rate is enhanced by coarse resolution. Because 
the ability of a given grid to support field growth through MRI stirring depends on the field 
strength itself when the resolution is marginal, excessive numerical reconnection can lead to 
a long-term suppression of the field intensity. 

This has implications for comparing physically significant quantities such as accretion 
rate or stress in simulations with differing magnetic field geometry. Without clear proof of 
convergence, the contrast between them as simulated with the same grid is as likely to be 
due to a contrast in effective resolution as it is to be due to a contrast in the real physical 
behavior of the different systems. 



5.5. Inflow equilibrium 



In shearing boxes one measures the properties of the MHD turbulence after many tens 
of orbits, when a quasi-steady state is clearly established. In global simulations one looks for 
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inflow equilibrium, a desirable property because it is the primary prerequisite for achieving a 
statistically stationary state. Any other state could, in principle, exhibit properties resulting 
more from transients due to arbitrary choices in the initial conditions than to the intrinsic 
physics of the system. Its meaning, however, requires precise specification in order to test 
how well a simulation matches this condition. The key question is how long to evolve a 
global simulation in order to reach a reasonable approximation to equilibrium. 

We have used several procedures for determining inflow equilibrium. One is to measure 
M(< R,t), the mass interior to radius i? as a function of time. A related procedure is to 
compute the characteristic inflow time from accretion rate and mass distribution, M/M, 
and the average inflow velocity (vr). This, in turn, can be compared to the inflow velocity 
derived using the time-averaged stress and steady state disk theory. Agreement indicates 
that the observed mass accretion rate is consistent with the observed angular momentum 
transport. For the fiducial run we find that all these measures give comparable results and 
show that the disk inside of -R ~ 30M is in inflow equilibrium. 



It is worth noting that iPenna et al.l (120101 ) proposed an alternative definition of inflow 



equilibrium time at radius R, specifically 2R/{v^ (rather than j dR/{vR))^ solely in terms 
of a general scaling argument. Because the magnitude of {vr) increases rapidly inward, there 
is a large difference between the local infall estimate R/vr and the actual integrated infall 
time. The Penna et al. value of the inflow equilibrium time is, for example, a factor of 7 
larger at i? = 20M than the actual inflow time J dR/ (vr) in the fiducial run. In fact, the 
inflow time estimated from 2R/ (vr) at the initial torus pressure maximum exceeds the time 
required to accrete the entire initial torus at the observed M. 

More generally, although an order of magnitude scaling argument may provide a rough 
rule-of-thumb for estimating the equilibration time, it does not provide a rigorous criterion. 
Only observed properties of the inflow can be used to test whether it is in equilibrium. After 
all, how a disk achieves a quasi-steady state is irrelevant; if the initial condition (miracu- 
lously!) agreed identically with the steady-state mass profile, it would be as good an equi- 
librium as one achieved after long computation. In fact, a commonly-used initial condition 
(one in which the magnetic field is dominated by large dipole loops; our 1-loop simulations 
are examples) behaves roughly in this way. Large turbulent magnetic stresses develop early 
in the simulation at the inner edge of the initial torus, pulling mass into the inner disk far 
more quickly than would happen at later, more representative, stages of evolution. Once 
the inner disk is filled with matter, however, it evolves toward an appropriate time-steady 
condition. The large transient magnetic stresses are simply a "jump-start" that allows a 
more rapid approach to the steady-state. 



As a final comment on this topic, we also wish to emphasize that achievement of in- 
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flow equilibrium does not, on its own, signal numerical convergence. In almost all of the 
global simulations we report here, the late-time state was one of approximate steady-state 
in terms of mass inflow within the inner disk, yet for none of them would we claim numerical 
convergence. 



5.6. Resolution and Numerical Dissipation Rate 

In simulations such as these, in which there is no explicit resistivity or shear viscosity, 
dissipation of turbulence occurs entirely due to gridscale discretization error (and any ar- 
tificial viscosity). This is true regardless of whether one uses an internal energy algorithm 
such as Zeus or an energy conservative scheme such as Athena. It follows that the numerical 
dissipation rate is affected by cell-size. Because short lengthscale dissipation is the primary 
means by which magnetic (and kinetic) energy is lost, the saturation strength of the tur- 
bulence depends directly on this dissipation rate. To understand how resolution affects the 
amplitude of the magnetic field therefore requires a closer consideration of gridscale losses. 



In the shearing box simulations of iHirose et al.l (120061 ) . the local rate of this sort of 
dissipation scales oc |J|^'^ for electric current density J. In other words, by Ampere's Law, 
numerical dissipation is roughly proportional to the cell- by-cell contrast in |B|. If the local 
scaling were exac% proportional to | J|, then the box- integrated dissipation would be likewise 
proportional to the box- integrated |J|. However, if there is any departure from linearity, it 
is necessary to be more careful in regard to the distribution function of current density over 
the box. 

The best data we have available for this purpose comes from the simulations of Shi et al. 
because their code explicitly tracked the several different kinds of gridscale dissipation. With 
remarkably little scatter, this data shows that the box-integrated magnetic dissipation rate 
Qtot oc Jtof, where Jtot is the box-integrated current density. Because this is numerical 
dissipation, not physical, it must depend on the cell-to-cell contrast in B, rather than on the 
actual current. Its most appropriate dimensional form is therefore 

n - ^ (^-^g"d)^^^ -1/4 /.ON 

Wtot - q — — p , [i-o) 

where g is a number of order unity. At is the time-step in the simulation, and p is the 
pressure. We evaluate (5-Bgrid) as (|V' x B|), where V denotes a curl operator defined with 
respect to cell-index rather than distance. As resolution improves, one would expect the 
cell-to-cell contrast to diminish, reducing Qtot- However, there are countervailing effects: in 
turbulence, new short lengthscale power can readily be generated by nonlinear transfer from 
longer wavelengths; and the time-step At also decreases. 
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In fact, the time-step At depends on spatial resolution in two ways: At = Ax/cg, and Cg 
depends on resolution because Shi et al. found that over their range of resolutions the total, 
i.e., gas plus radiation, pressure had not yet converged; it scaled roughly oc (Ax)^, with k ~ 
—0.45. At fixed surface density, hydrostatic balance implies that Cg oc p, so At oc {AxY^'^. 
Thus, we expect that the overall resolution scaling of Qtot should be oc {5Bgnd)^^'^{Ax)~^^^ . 
Comparing the Qtot versus (5-Bgrid) relations found in the Shi et al. simulations, we find 
that the actual exponent of Ax is ~ —1.3, in very good agreement with the dimensional 
argument leading to Equation (ITSI) . 

When the magnetic field achieves its saturation strength, Qtot balances the fluctua- 
tion power delivered to the grid scale by nonlinear coupling with longer lengthscales in the 
MHD turbulence. At resolutions currently feasible, these couplings are non-local in the sense 
that the shortest wavelengths couple significantly to essentially all longer lengthscales; how- 
ever, there is reason to think that in genuine MHD turbulence, where the dynamic range 
between stirring scale and dissipa tion scale is far greater, the dissipation scale actually de- 



couples (ILesur fc Longarettil 120101 ). Because all our simulations are in the regime in which 
all lengthscales couple directly to the grid scale, we write the rate at which longer scales 
couple to the dissipation scale as 

i^nonlin = sQ{6BI,) /S'^ , (19) 

where {6B^^^) is the total variance in the magnetic field. Nonlinear energy transfer occurs 
on the eddy turnover, i.e., dynamical, timescale and involves the entire fluctuating power. 
Equating i^noniin to Qtot yields the prediction 



SB. 



tot 



in 



-p m+3/4 f5/2 

s Sfi2Axo(Ax/Axo)"'+i-3*^/4 



l/(2m-l/2) 

(20) 



Here we have introduced the adiabatic index F, column density S, orbital frequency fl, char- 
acteristic pressure po, and grid resolution Axq for which p = pq. We have also introduced the 
quantity /grid, the ratio between the characteristic amplitude of magnetic field fluctuations 
on the grid scale and on the stirring scale. 

As the cell-size Ax decreases, we can expect /grid to decrease. However, only when 
convergence is reached must it scale oc Consequently, convergence in the mag- 

netic field intensity may not be monotonic. Using the data in our stratified shearing box 
calculations, we can measure how far toward convergence we have progressed in this re- 
gard by computing /grid- Combining the Simon et al. and Davis et al. data, we find that 
/grid falls from 0.42 for the simulation with 8 cells per scale height to 0.26 when there are 
128 cells per H. In other words, for these simulations, in which k = and the resolution 
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improves by a factor of 16, /grid falls by only 40% so that, if the index m = 0, the ratio 
/g£/(Ax/Axo)'"+^-''=/^ would rise by almost a factor of 5 and the saturated field strength 
would fall by a factor of 25. That the saturated magnetic field in fact increases, but rather 
slowly, over this range of resolutions suggests that — 1 ^ "^i < —0.5; in other words, the 
nonlinear energy transfer rate increases as the gas pressure rises relative to the magnetic 
pressure. With this dependence of the nonlinear energy transfer rate on the plasma /3, the 
saturated field strength increases as the gridscale fluctuations become smaller relative to the 
total variance. At still flner resolution, the scaling with Ax of both /grid and the nonlinear 
coupling could change, but at least in this range of resolutions for these algorithms, the fleld 
intensity can be expected to grow with further reflnement in spatial grids. 

Thus, the action of gridscale dissipation creates a direct relation between resolution 
and the nonlinear processes that determine the magnetic fleld intensity and Maxwell stress. 
Because improving resolution also increases short lengthscale fluctuation power, the level of 
gridscale fluctuations diminishes only slowly with improving resolution. If, as appears to 
be the case in the simulations we have studied, nonlinear coupling between large scales and 
small grows modestly with an increasing ratio of gas pressure to magnetic pressure, the net 
result is a slow increase in saturated magnetic intensity with improving resolution. Processes 
like these account, at least in part, for the more stringent requirements for resolution that 
we have uncovered. 



5.7. DifRculties of convergence testing 

Although convergence testing is an accepted approach for gauging the reliability of 
numerical simulations, there are several issues that make it extremely difficult to apply 
to global accretion simulations. Some of these difficulties are practical; others are more 
fundamental. 

One of the fundamental difficulties arises because accretion is driven by MHD turbu- 
lence, and the interaction of dynamics at widely-differing lengthscales is intrinsic to the 
nature of turbulence. Large dynamic range in spatial scales generically demands very high 
resolution, of course. The difficulties associated with achieving convergence are especially 
compo unded for MRI-driven MHD turbulence for which stirring occurs over a wide range of 



scales. iLesur fc Longarettil (120101 ) argue that a clearly defined inertial range separating the 
stirring scale from the dissipation scale still has not emerged in the shearing box, even at 
the highest resolution available to their spectral code. Because the total Maxwell stress is 
dominated by field energy on the largest scales, one may hope that its magnitude becomes 
reasonably well-determined even while the resolution remains insufficient to describe details 
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of the short-wavelength dissipation, but that remains to be demonstrated. The previous 
subsection discussed some of the problems to overcome in order to achieve that decoupling. 

Among the practical considerations is that MHD turbulence-driven global accretion 
simulations must be three dimensional. Achieving even modest resolution in three dimensions 
requires a large number of cells. Whenever the science pushes computational capacity to its 
limit, it is infeasible to carry out additional simulations at higher resolution than the level at 
which one begins to see interesting effects. For this reason, attempts are sometimes made to 
"test convergence" by contrasting one set of results with those found at coarser resolution. 
Although this may be the only available procedure for simulations carried out at the largest 
feasible resolution, it is an approach that must be used with caution. If the gridscale at every 
resolution tested remains too coarse to be sensitive to a physical mechanism, the results will 
give the appearance of convergence, even while missing the very existence of a mechanism 
that may be very important. This motivates this paper: using the better-resolved shearing 
box simulations to develop diagnostics that allow a global simulation to be evaluated directly. 

Finally, there is the difficulty of measuring the effects of truncation error and the rate 
at which it is being reduced for simulations carried out at relatively coarse, but practical, 
numerical resolutions. Strictly speaking, every algorithm has its own properties, particularly 
with respect to truncation error. Truncation error is characterized by the amplitude of the 
error, the effective convergence rate at a given resolution, and the asymptotic convergence 
rate. The latter is a fundamental property of the method, but in practice one almost never 
reaches the asymptotic regime, except in select test problems or for one- dimensional prob- 
lems where exceptionally fine resolution may be achieved. At resolution levels reached in 
typical three-dimensional simulations, the observed convergence rate need not agree with the 
asymptotic convergence rate of the algorithm. Also, the convergence rate alone provides no 
information about the amplitude of the truncation error at a given resolution. Again, by de- 
termining properties that are characteristic of well-developed MRl-driven MHD turbulence 
(e.g., {Bj^/B'^)), we sidestep this difficulty. 

Both global and local simulations have made use of several different algorithms (finite- 
difference: Zeus, GRMHD; finite- volume conservative: Athena, HARM3D), as well as varia- 
tions with in the same general algorithm (e.g., different flux solvers). The relative merits of 
different algorithms will be reflected in the value of the Q parameters required to produce 
good results. For the shearing box results in § |3l we found that the use of either Zeus or 
Athena appears to make relatively little difference; our diagnostics of convergence depend 
upon resolution in very similar ways. In § E] we compute the number of grid zones that 
might be required for a future well-resolved global simulation. These numbers can vary for 
one algorithm compared to another, according to the Q values that prove necessary. 
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5.8. Spatial resolution requirements and published global simulations 



There is no single number by which spatial resolution may be defined. Even to specify 
the raw cell count requires three numbers, one for each dimension, and that is not sufficient 
because cell size depends on the extent of the simulation volume. Moreover, what really 
matters is not so much the absolute size of the cells as their size relative to the natural 
lengthscales of the problem. For global simulations, the radius R is the primary scale relevant 
to the radial direction, so the appropriate measure of cell size is AR/R. In the azimuthal 
direction, there are two natural scales: radians and 27rv A(I,/^- The former is due to the 
necessity of resolving a wide-range of non-axisymmetric modes in the turbulence; the latter 
has to do with adequate resolution for the fastest-growing linear non-axisymmetric modes. In 
the vertical direction, there are likewise two scales: the density (or pressure) scale height, H, 
and 27tvaz/^- Quantities based on the Alfven speed always require time-averaging because 
the field evolves rapidly over the course of the simulation; depending on the treatment 
of thermodynamics in the simulation, the vertical scale height may also require averaging. 
Therefore, to align resolution metrics against simulations in the literature requires translation 
of all cell-size specifications into these five natural units. 

We have done this for a selection of published simulations in Table |H In addition to cell 
size data, we list the radius of the pressure maximum in the initial torus used, and the simula- 
tion length in terms of orbits at that radius. The data from the various simulations reported 
here with the L, M, and R grids are also given in this table. The add itional simulations 



used are: the high-resolution poloidal-field simulation of iHawlev &: KrolikI (|2002[). w hich we 



De Villiers et al 



(j2003|), called 



label HK02; the high-resolution Schwarzschil d simulation of I . 

KDOb; the quadrupolar field simulatio n from Beckwith et al] ~( 2008b ). designated QDO; the 
high- resolution thinnest simulation of Noble et al. ( 201o[ ). which they labeled ThinHR; the 



Schwarzschild simulation reported by IShafee et al.l (l2008il. identi fied as "Shafee" in the ta- 
ble; and the fiducial simulation A0IIR07 from lPenna et al.l ( l2010l ). here called "Penna". For 
many of these, the published paper did not report full information. For those for which we 
have access to the data, we have computed what we needed; for the remainder, we have filled 
in only those entries derivable from the publications. Where this required reading values off 
graphs, we have added a question mark. All but one of the simulations listed in the table 
used an azimuthal range of 7r/2; the Shafee simulation used a wedge only ir/A in angle. In 
these other simulations, unlike the Zeus simulations described in this paper, the radial grid 
has a fixed AR/R. They also differ in that, unlike our new Zeus simulations, which used 
cylindrical coordinates, they used spherical coordinates. Because Az = RAO in spherical 
coordinates, the radial gradients in {Q^) seen in the Zeus simulations are greatly reduced. 



As the data in this table demonstrate, most global simulations performed hitherto have 
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had effective resolutions rougfily at tfie level of our M grid with the standard azimuthal 
resolution, while the Penna simulation most closely resembled our twoloop-lOOO-mlp simu- 
lation, which used our poloidal M grid but coarser azimuthal resolution. As a fraction of the 
local radius, the radial cell sizes have been generally in the range ^ 0.01-0.025, similar to the 
0.017-0.055 range of our M grid; the only exceptions in this list are Shafee {AR/R = 0.0065) 
and ThinHR {AR/R = 0.004), both of whose radial resolution was closer to our R grid. In 
terms of the scale height, the vertical cells of most of these simulations are in the range 0.02- 
0.14, again matching our M grid; by this measure, Shafee is similar to all the others, while 
ThinHR is the only one substantially better, with Az/H — 0.0086. These values match the 
range explored in the stratified shearing box simulations, where the coarsest we discussed 
had a vertical cell 0.125i7 thick, and the thickness of the finest vertical cell was O.OOTSiJ. In 
radian measure, all of these simulations except Penna used a cell with A0 = 0.025, identical 
to that in our standard M grid; Penna's azimuthal cell was twice as large. Our ability to 
be quantitative about (Qz) and (Q^) is limited: we kept insufficient data from HK02 to 
compute these diagnostics, and we lack access to data from Shafee and Penna. However, on 
the basis of what the published papers say, it is clear that {Qz) in most of the global disk 
literature has been ~ 5-8, while (Q^) has been ~ 20. By these measures, at ~ 25, ThinHR 
was a lone standout in terms of (Qz), but in the same range as all the others in terms of 

{Q<t>)- 

Because our M and R grid simulations at best went only part way toward convergence 
as measured by any of our diagnostics, we expect that all these simulations likewise fell 
short of convergence. For those cases in which the initial magnetic field was in a single 
dipole loop (HK02, KDOb), the {Qz) values suggest that the {B\/B'^) measure is in the 
neighborhood of ~ 1/2 its saturated value. For those in which the initial magnetic field 
was in two loops (QDOb, Shafee), the shift we have seen in 2- loop simulations relative to 
1-loop suggests that {B\/ BV) is somewhat smaller. The Shafee simulation also has large 



Table 4. Cell-sizes 



Name 


AR/R 


A<t> 




Az/H 




Rmax 


Orbits 


L Grid 


0.033-0.111 


0.025 


7-10 


0.0333-0.1 


1-3 


35 


33 


M Grid 


0.017-0.055 


0.025 


~ 10 


0.016-0.08 


2-10 


35 


54 


R Grid 


0.004-0.0067 


0.025 


6-12 


0.02-0.1 


3-10 


35 


28 


HK02 


0.008-0.015 


0.025 




0.02-0.07 




20 


15 


KDOb 


0.024 


0.025 


19 


0.06-0.14 


6 


25 


14 


QDO 


0.024 


0.025 


26 


0.08 


8 


25 


14 


ThinHR 


0.004 


0.025 


18 


0.0086 


25 


35 


12 


Shafee 


0.0065 


0.025 




0.03 




35 


8 


Penna 


0.013 


0.05 




0.064 


~5? 


35 


22 
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aspect ratio cells (although Ai?/ {RAcp) is similar in ThinHR, the latter's much finer vertical 
resolution appears to compensate). In Penna, the combination of a four- loop configuration 
and coarser azimuthal resolution seems particularly challenging. In our simulations such 
conditions cause {B\/ B"^) to fall by a factor of 2 or more. Such a reduction in the relative 
magnitude of the radial component of the magnetic field would indicate that the turbulence 
is not fully developed, that the field intensity is rather lower than the converged value, and 
that the stress is smaller than in saturation because it is proportional to Br. 



6. Conclusions 

Although global disk simulations have explored many important aspects of the accretion 
process, their quantitative reliability remains uncertain. In this paper we have made use of 
high resolution local shearing box simulations to develop four diagnostics by which one 
may gauge how closely a given simulation approaches to fully developed MRI-driven MHD 
turbulence. We have then examined how those diagnostics carry over to global simulations. 
Establishing this connection is an important step in relating local simulations to global. Local 
simulations will always be able to include more physics and use better effective resolution 
than global, and by means of this mutual calibration, information from shearing boxes can 
help guide and interpret global models. 

These four are: {Qz)-, the number of vertical cells across a wavelength of the fastest- 
growing poloidal field MRI mode; {Qy) (or {Q^) in global simulations), the number of az- 
imuthal cells across a wavelength of the fastest-growing toroidal field MRI mode; {B'^/B'^), 
the ratio of energy in the radial magnetic field component to the toroidal component; and 
amag, the ratio of the Maxwell stress to the magnetic pressure. Only the first has seen sig- 
nificant use previously, and we extend its utility to gauge spatial resolution fo r nonlinear 



behav ior of the MHD turbulence as well as linear growth of the MRI. Whereas ISano et al. 



(120041 ) found that a minimum (Q^) of 6-8 is required in order to describe poloidal MRI linear 
growth, we find that the prerequisite for simulating nonlinear behavior is more stringent and 
couples poloidal resolution to azimuthal. If the analogous azimuthal diagnostic {Q^) > 20, 
then {Qz) > 10 is necessary; if {Q^) is smaller than that, still larger values of {Qz) are 
required. 

Using these diagnostics, we find that all the global simulations done to date are under- 
resolved, particul arly in the (p direction. Only one simulation in the literature (ThinHR, 



Noble et al.ll2010l ) comes close to adequate poloidal resolution, but even it does not meet the 
azimuthal standard. Our tests varying the azimuthal cell-size showed that it can be impor- 
tant to the development of the poloidal, as well as the toroidal, magnetic field. Achieving 
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adequate azimuthal resolution is made especially difficult when the disk is thin for two rea- 
sons. Thinner disks require smaller vertical cell dimensions, but avoiding the deleterious 
effects of large cell aspect ratios then demands cells still smaller in the azimuthal direction. 
In addition, va tends to diminish with the smaller sound speeds seen in thinner disks; smaller 
A0 is then necessary in order to achieve an adequate value of Q^. 

Regarding the initial conditions, we have found that the additional short lengthscale 
fluctuation power that is a concomitant of more complicated magnetic field geometry places 
stronger demands on spatial resolution. The additional magnetic reconnection associated 
with attempting to describe a more complicated field geometry on a given grid can weaken the 
field sufficiently that MRI stirring is curtailed and the field intensity remains artificially low. 
Until global simulations are adequately resolved it will be difficult to distinguish numerical 
from physical effects arising from different initial field geometries. 

We can make an estimate of the resolution required for a global torus evolution using 
Equations (fTTj) and (fT3|) . To avoid the problem of decreasing H with decreasing R that is 
found in cylindrical coordinates, we assume a spherical-polar grid (r, 9, 0). To keep Ar = rA0 
we use logarithmic spacing in r. We also assume the radial grid spans a factor of 100 from 
the inner to the outer boundary, and that the (f) grid spans only 7r/2. We next assume 
/3 = 10, 13 J 13 = 50, H/R = 0.1. Our target values are {Q^) = 10 and (Q^) = 25 (one's 
preferred target Q values may well be algorithm dependent). With these assumptions, the 9 
grid, if equally spaced, must have ~ 900 zones if equally spaced from to vr. In a practical 
simulation, however, one could increase A9 away from the midplane, reducing the number 
of zones required, perhaps by as much as 50% (450 zones). The grid requires 200 zones 
and the r grid 600 zones. For H/R = 0.1, the number of cells (600 x 200 x 450) is not 
too much greater than what has already been used. Of course, reducing H/R would further 
increase the number of required zones proportionally in each of the three dimensions. Higher 
resolution also reduces the Courant-limited timestep, making it challenging to evolve the disk 
for a large number of orbits. In this estimate we assumed a smooth distribution of field; 
intermittency in the field distribution may change the required resolution from this estimate. 
Typically there is significant spatial variation in (3 so that a density weighted Q value can be 
larger than a volume-weighted value. On the other hand, complex initial field distributions 
may require more zones for adequate representation. Future simulations done at resolutions 
approaching this estimate can test these ideas. 

Although we find that certain values of the resolution parameters indicate well-resolved 
turbulence, lower values of those parameters don't mean that turbulence necessarily decays 
rapidly; as we have described, it can persist for long periods of time with lower amplitude 
and reduced stress. Even under-resolved global simulations can have nonzero stress levels 
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and accretion rates. A seemingly long duration of accretion is therefore not a guarantor of 
convergence. Qualitative conclusions can usefully be — and have been — obtained from such 
simulations, but quantitative stress levels are likely to be undervalued. 

Lastly, we have explored the question of the meaning of "inflow equilibrium" in the 
context of global simulations. Analysis of several diagnostics can determine the existence 
and extent of the region of the accretion disk that is in inflow equilibrium. One can show 
that there are at most weak long-term trends in the radial mass distribution, that the 
time-averaged accretion rate is relatively constant as a function of /?, and that the observed 
accretion rate is consistent with the observed angular momentum transport as computed from 
the steady-state disk equation. Even when the disk is in time-averaged inflow equihbrium, 
there can be large short-term fluctuations in the mass distribution. Moreover, whenever 
there is only a finite mass reservoir in the initial condition, a sufficiently long simulation 
cannot, by definition, support even a statistical steady-state for longer than the inflow time 
from the initial half-mass radius. 

In summary, this work provides guidance for future global simulations, both in terms 
of resolution and evolution time, to approach what is needed for quantitative conclusions 
about accretion disk dynamics and structure. 
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Fig. 1. — Initial torus and field configuration for the two-loop simulations. The outermost 
contour line marks the boundary of the initial torus. The remaining contours are the mag- 
netic field lines. 
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Fig. 2. — Ratio of «, the volume- averaged Maxwell and Reynolds stress to volume- averaged 
pressure in a set of four stratified shearing box simulations. The simulations use 8 (lowest 
solid line), 16 (dotted line), 32 (dashed Une), and 64 zones (solid line) per scale-height H. 
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Fig. 3. — Distribution function for the time- and volume- averaged ratio of the Maxwell stress 
to the magnetic pressure, a„iag taken from a set of four stratified shearing box simulations 
using 8 (left-most solid line), 16 (dot-dashed line), 32 (dashed line), and 64 zones (solid line) 
per scale-height H. As resolution increases the peak of the distribution of a^ag shifts to 
higher values. 
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Fig. 4. — Ratio of volume-averaged radial to toroidal field energies, {Bl)/{By) for stratified 
shearing box simulations using 8 (solid line), 16 (dotted line), 32 (dashed hne), and 64 zones 
(sohd line) per scale-height H. This ratio shows a systematic increase with resolution. 
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Fig. 5. — Ratio of the volume- averaged radial field energy to toroidal field energy, (Bl/By) 
as a function of resolution in stratified shearing box simulations. Resolution is measured in 
number of grid cells per scale-height H. 
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Fig. 6. — Evolution of the three components of the magnetic energy as a function of time in 
the fiducial simulation, twoloop-lOOO-mr. 
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Fig. 7. — The mass interior to three radii (i? = lOM: soUd, R = 15M: dotted, R = 20M: 
dashed) as a function of time in the fiducial simulation. Mass units are fraction of the total 
initial mass. The time averaged accretion rate, M, is 1% of the initial mass every 3700M in 
time. 
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Fig. 8. — The inflow velocity for the fiducial run derived from the vertically and azimuthally 
integrated accretion rate and density (solid line), along with the infall velocity derived from 
the vertically and azimuthally averaged Maxwell stress, pressure and accretion rate using a 
steady state disk approximation (dashed line) with j„ = jisco- The dot-dashed line is the 
steady state inflow velocity derived using = 0.985jisco- The data are time-averaged from 
i = 4-5 X lO^M. 
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Fig. 9. — Evolution of the magnetic energy in the fiducial run (solid line) compared to the 
magnetic energy in a stratified shearing box simulation with 16 zones per H (dashed line). 
Time is orbits at the initial torus pressure maximum or in terms of the shearing box and 
the magnetic energies are normalized to their peak value. Both models show a period of 
rapid field amplification followed by a peak and a decline to a longer-term value that slowly 
declines. The magnetic energy increases after 50 orbits in the shearing box, but not in the 
global model. 
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Fig. 10. — Radial dependence of the principal diagnostics for the 2-loop simulations with 
the standard azimuthal resolution. The curves are labeled with model name. The data were 
averaged over azimuth and height with a density weighting and then averaged in time from 
2-3.4 X lO^M for all three simulations. (Top) (Middle) (Q^), (Bottom) (Q^). 
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Fig. 11. — Radial dependence of the principal diagnostics for the 2-loop simulations with 
varying azimuthal resolution. The lines are labeled with model name. The data were av- 
eraged over azimuth and height with a density weighting and then averaged in time from 
1.8-2.5 X lO^M for all three simulations. (Top) (Middle) (Q^), (Bottom) (Q^). 
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Fig. 12. — Evolution of three simulations using different numbers of zones to span the (f) 
domain: 32, 64 (dashed line) and 128. The initial conditions were the same: two magnetic 
field loops with average strength P — 1000. 
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Fig. 13. — Radial dependence of the principal diagnostics for the 1-loop simulations with 
the standard azimuthal resolution. The curves are labeled with the corresponding model 
name. The data were averaged over azimuth and height with a density weighting and then 
averaged in time from 2-3.8 x lO^M for all three simulations. (Top) (Middle) 
{Q,). (Bottom) (Q^). 
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Fig. 14. — Correlation between (Q^) and {Bj^/B'p. The points shown by filled circles in this 
figure are drawn from the same set of times in three different simulations: twoloop-betalOOO- 
Ir, twoloop-betalOOO-mr (the fiducial global simulation), and twoloop-betalOOO-mhp. All are 
density-weighted averages, but they arc taken at different radii within the inner disk. The 
points shown by the other symbols are time-averages of data from the midplane region of 
the shearing box simulations. The different symbols within this category denote different 
ranges of (Qy): 12 < (Qy) < 25 (x's); 25 < (Qy) < 50 (squares); (Qy) = 98 (diamond). The 
dashed line has a slope of 0.65. 
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Fig. 15. — Time-averaged and shell-integrated Maxwell stress profiles for the three 2-loop 
simulations differing only in poloidal grid (Top) and the three 2-loop simulations differing 
only in azimuthal grid (Bottom). Curves are labeled by their model name. 



